17/05/06 16:39:10.74 NK3Wh0hW.net
当直スレに以前書いたクイズ
A君の彼女は女子大生、B君の彼女は女子高生。
A君は女子大生の方が胸が大きいと主張して次のようなデータを出した。
女子大生100人と女子高生100人の胸囲を測定して
前者が平均82 , 標準偏差3
後者が平均81 , 標準偏差3
であった。
t検定でp = 0.0194と有意差があったとA君は勝ち誇っている。
さて、女子高生の彼女をもつB君は平均1cmの差では実地臨床w上、意味がない、
5cm以上の差があってこそ体感上有意と言えると主張。
問.
B君の女子高生彼女がA君の女子大生彼女より5cm以上の巨乳である確率とその95%信頼区間を述べよ。
新薬投与群100人でHbA1c7.2(sd=0.3)、従来薬100人でHbA1c7.1(sd=0.3)で統計的に有意差ありという医学ネタにすると面白くないね。
シミュレーションでの計算例は
ｽﾚﾘﾝｸ(hosp板:756-757番)
23:卵の名無しさん
17/05/06 17:10:55.66 NK3Wh0hW.net
MCMCでシミュレーションしてみた
Inference for Stan model: boobs.
4 chains, each with iter=15000; warmup=5000; thin=1;
post-warmup draws per chain=10000, total post-warmup draws=40000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
d 0.990 0.024 4.251 -7.352 0.974 9.308 31862 1
JD 0.175 0.002 0.380 0.000 0.000 1.000 28355 1
JK 0.078 0.002 0.268 0.000 0.000 1.000 26103 1
dは胸囲の差
分布をグラフにすると
URLﾘﾝｸ(i.imgur.com)
JDは女子大生のバストが女子高生より5cm以上大きい確率、
JKは女子高生のバストが女子大生より5cm以上大きい確率。
新薬の切り替えたときに臨床的有益性がどれくらいあるかの確率が計算できる。
24:卵の名無しさん
17/05/06 18:42:59.31 gFsvkygs.net
>>22
今日の放送大学の講義で、pって何だったんですかって話になって、
「元の確率分布から非常に離れた値を取る確率」
という説明で、
「元の確率分布が正しい確率が1-pになるわけじゃないですよ」
と説明されていた。
有意差検定とベイズ統計学の差が出ているのかなと思いました。
25:卵の名無しさん
17/05/06 21:48:23.39 ++fNBVM2.net
>>24
MotuluskyのIntuitive Biostatisticsに
ｐ値がでたら、帰無仮説は何かを考えよ という趣旨のことが書いてあった。
リスク比RR＝１ （またはリスク差RD＝０）を帰無仮説にしてｐ値を算出するけれど
RR＝１に帰無仮説を設定する必然性はない、例えばRR＝１．５に設定してｐ値が算出できるという記載が
Rothmanの疫学(introductionの方)にあって
そのグラフｐ値関数を再現できるまでずいぶんと悩んだ（寝当直の暇つぶしになった）。
26:卵の名無しさん
17/05/06 23:44:38.35 ++fNBVM2.net
>>24
その誤解がアニメになっている。
URLﾘﾝｸ(www.youtube.com)
p.value functionを検索していたら
URLﾘﾝｸ(sphweb.bumc.bu.edu)
に掲載されていた。
27:卵の名無しさん
17/05/07 00:46:48.66 1tZb4eYU.net
>>26
かわいいアニメですが、日本語のアニメと同時には聞き取れませんね。今日は無理か、
ガンダム観たいし（笑）、明日にします。
後、やっぱり、有意性検定で書かれた論文が多いので、便利なベイズ統計と二足の
ワラジになりそうですね。
まだ、リスク比とハザード比の差が分からないのですが...。
28:卵の名無しさん
17/05/07 01:13:01.15 glcxZk26.net
>>27
リスクは一定、
ハザードは時間の経過と共に変化する。
外科治療と内科治療でのハザード比は一定というモデルは
周術期と安定期では一定でないと思う。
リスク比とハザード比をひっくるめてRate Ratioと呼んだりする。
RRはどちらの略か区別がつかないこともしばしば。
比例ハザードモデル=ハザード比一定
ならリスク比と区別がつかなくてもいいんだろうと理解している。
29:卵の名無しさん
17/05/07 15:19:34.12 1tZb4eYU.net
>>28
つまらない、ピエロみたいな話にレスしてくれてありがとうございます。
コホート研究と、一回の2x2の検定の差は分かります。
で、コホート研究の確率分布関数が指数関数という仮定に基づいて
いませんか？
それは、統計学者にとっては致命的ではないのでしょうか？
私は、それががん治療における
「最終命題」
じゃないのかと思うわけです。
「つらい」の変換で「痛い」がでないんですか？
「つらい」は日本人的には「痛い」なんですよ。
先生が英語で話したい理由は分かります。世界中に関する影響力が
違いますよね。
英語で聞き、英語で表現する、それが確率的に影響力最大なんですよね？
「日本人は苦労して、世界に届いている」
どうしてそれだけの能力があるのに、世界で認められないのか
そういう話ですよね。
30:卵の名無しさん
17/05/07 15:22:50.66 1tZb4eYU.net
確実に臨床統計は、
「日本人の生死とQOLを握っている」
わけでしょ。
日本人が複数である以上、それを否定する理論があれば、間違いですよ。
31:卵の名無しさん
17/05/07 15:26:28.75 1tZb4eYU.net
医師が努力すれば、
「医師法第一条に定められた、公衆衛生の利益が守られるわけじゃない」
という時代に、
利益最大の医師会や病院団体と、
厚労省の利益は合わないでしょ？
32:卵の名無しさん
17/05/08 14:46:04.38 vkqJQ/ib.net
>>30-31
ごめんなさい、
他の板の関係で、数学板やら物理板の関係でこんなことを書きましたが、
学問は中立ですよね。中立ゆえの権利があるはずですね。
と思っているわけです。
数学板でモンティ・ホール問題に疑問を呈して、恣意的な確率操作に対して
意を表明していますが、
人間は、理論より納得でしょ。それを得られるかどうか、やってみているわけです。
ダメなら、呪術をmedicineと呼んでいるアフリカ諸国と同じですよ。
33:卵の名無しさん
17/05/08 14:52:32.64 vkqJQ/ib.net
URLﾘﾝｸ(detail.chiebukuro.yahoo.co.jp)
日本語で、「ハザード比、オッズ比、リスク比」を検索すると
こういうyahoo知恵袋が一番上に出てくるわけで、
情けない
と思うわけです。
誰が日本の臨床統計に対して責任を持っているのか、そういう学会があれば、その検索結果が
トップにでるように努力すべきなのではないですかね？
ネット時代に対応していない学会は捨てられるか、利益団体に利用されるだけかな？
34:卵の名無しさん
17/05/09 15:02:02.15 pn7nUozD.net
東京馬鹿大で、
「医学生が医師を刺した」
という事象はどう説明するのかね？
明らかに、これは確率論・統計学的な話と言えないかな？
まともな臨床統計を目指すならコメントすべきじゃないし、
言いたいことがある人なら、コメントすべきじゃないかな。
日本語は使いにくいかな。
35:卵の名無しさん
17/05/09 15:05:27.95 pn7nUozD.net
もっと言えば、
「日本人であること、日本国籍を持っていることは有利なのか？」
って話だと思うよ。
反論したければ、公的保険に対する期待値を示すべきだね。
36:卵の名無しさん
17/05/13 14:05:24.66 Z5ocLfe6.net
AI=機械学習が人間の決定論に対して、鏡になってくれているわけかな。
EBMで、医療すら確率論的空間になってしまったが、
「人間という機械は、するかしないかの二択でしかあり得ない。」
だね。
「確率を提示するのが医師」
なのかな？ 論理学ならすべての事象を集合論で網羅して、真偽を与える
わけなんだけれど。
私はおかしげな当時の決定論に疑問を持って、大学院は拒否したわけだけれど、
多くの患者さんの犠牲の下に、私は地位を確立したかったのか、それを拒否したのか、
その感覚的判断は正しかったのか、バカみたいな論文を書き続けて今の日本の
医学的な指導者になっているのか、結局、
「屍の上に学問は成り立っているわけ」
でしかない。それをやった奴は手を血に染めているわけだ...。
その実験をするためには、それを誘導するための意志があったわけだろ。
二重の事後確率を決定するわけだよね。
37:卵の名無しさん
17/05/13 14:13:03.97 Z5ocLfe6.net
だいたい、普通の臨床医は
「即断即決」
を求められているから。
EBMは嫌いなんじゃないかな、と思う、私もそうだから。
外来の決められた時間で
「決定論的結論」
を求めている患者さんには、この考察にかかる時間は
「不利益」
なんだよね。
確率論は熱力学と同じで、時間的問題を無視している。
サイコロはいつ投げられても同じ、
「コンビニの店長は、時間という価値と確率論的価値
はどう判断しているのか？」
簡単なORにおける在庫管理じゃないよね。
38:卵の名無しさん
17/05/14 13:05:43.27 75s6edjB.net
>>37
治療による診断はさしずめ、シュレーディンガーの猫かな？
39:卵の名無しさん
17/05/14 13:10:25.01 75s6edjB.net
>>37
天気予報にも降水確率つきで発表されるから、正診確率とともに診断を伝えてもいいんじゃないか？
事前確率は往々に経験則もしくは主観的設定だけど。
40:卵の名無しさん
17/05/14 16:44:04.43 a+Ds8MZq.net
スレチと思うが、AIは医者業務（手技や承認作業を除く）の何を置き換える事ができるかな？
という問いに対して医者から「絶対これは譲りたくない」と心理的抵抗が強い業務こそ、AIを開発することによって医者減らしに貢献すると思うがどうかな？（救命救急士が行う気管内挿管のように、医者の既得権益？に踏み込む事だろうから。）
41:卵の名無しさん
17/05/14 17:56:55.21 AzIdUf7d.net
>>40
画像診断（ミクロの画像診断、病理も含む）だな。
42:卵の名無しさん
17/05/15 00:11:14.10 j+qy95C6.net
あの話題の線虫によるがん検診の話がどのぐらいの信憑性があるのか、
知りたいなぁ。
43:卵の名無しさん
17/05/15 00:24:43.07 j+qy95C6.net
>>40
AIの能力次第だと思うんだけれど...。
大量のデータから診断や治療法を探す検索系と、
機械学習で得られた成果を組み込んだ、画像判断ソフト
みたいな判断系と、熟練の技術者のデータを入れた
自動車組み立てロボットみたいな手技系が
できるんじゃないのかな？
44:卵の名無しさん
17/05/20 14:35:39.61 1Y/5oN2J.net
土曜の午後は放送大学の無料放送で
「ベイズ統計学の復習」
の予定だったのですが、数週間できませんでした。
今日は見ます。
45:卵の名無しさん
17/05/20 14:48:47.08 1Y/5oN2J.net
>>44の続き
今日も無慈悲に
飛ばしている講師の先生。
不偏分散の話もどこかの回で無視していたのだが、
nで割るのとn-1で割るのと
たいした違いがない。分かるんだけれどね...。
URLﾘﾝｸ(eman-physics.net)
にきれいな証明があるからご参考までに。
46:卵の名無しさん
17/05/20 14:58:43.46 1Y/5oN2J.net
R言語で積み重ねられた、
統計学に必要なパッケージと
ランダムに作る、100万以上の仮想的なデータ
があれば、
統計学の先人達も、そりゃ対数表に頼っていたニュートンやら
昔の数学・物理学者も1,000円ぐらいの電卓に負けるかもね。
考慮すべき技術的な基盤技術が違う...、
事前確率分布の関数ことが人間の決められる時代なのかな。
47:卵の名無しさん
17/05/20 15:13:07.04 1Y/5oN2J.net
あれ、今日もp値の評価の仕方かな。
最後に有意性検定の否定かな？
まず、帰無仮説の否定...
まあ、母集団とサンプル群の値が同じだとは言えないよな。
サンプルを取る前から偽...
そりゃあ、恣意的なのかもな。
Jargon!!!!
笑います。
ずごい！ 刑事事件になってきたぞ。
アリバイが検定に導入されている。
アリバイがないことで有罪の確率が変わるのはおかしい、
まあそういうことかな。
48:卵の名無しさん
17/05/20 23:59:15.59 ZMCZkpyW.net
>>47
親子鑑定にベイズが使われたりするんだよな。
理由不十分の原則とかに陪審員が騙されるんだよね。
49:卵の名無しさん
17/05/22 18:02:30.07 bGYMjDtG.net
npoでお金悩み相談。
日々の生活での返済、お支払いでお悩みの方。
急な出費などで、今月の生活費が足りない方。
多重債務、ヤミ金、家賃滞納でお困りの方。
お金に関するお困り事や法的トラブル等HPに記載以外の事でも、お気軽にご相談下さい。
東京、神奈川、千葉、埼玉にお住まいの方は優遇です。
詳しくはHPをご覧下さい。
ｴｽﾃｨｰｴｰで検索
50:卵の名無しさん
17/05/24 19:20:21.54 olQQP7cX.net
女医が担当医の方が生存率が高いという論文
（１）Comparison of Hospital Mortality and Readmission Rates for Medicare Patients Treated by Male vs Female Physicians
URLﾘﾝｸ(scholar.harvard.edu)
若い医師が担当医の方が生存率が高いという論文
（２）Physician age and outcomes in elderly patients in hospital in the US:observational study
URLﾘﾝｸ(www.bmj.com)
この二つの論文から
（３）若い女医が担当医のときが生存率が最も高い
（４）老いた男性医師が担当医のときが生存率が最も低い
という結論を導きだせるか。
51:卵の名無しさん
17/06/10 16:09:52.77 6oyq1OSX.net
ありゃ、しばらくアクセスしていないうちにスレ落ちちゃったわけで
申し訳ないです。
確率論におけるヒルベルト空間での議論で
はあ、有意差検定というのは、
10,000人のデータを真にしたとき
の30例のデータが優位であるかないかの話で、
10,030例の
検証に至っていないわけです。
臨床統計に至る100例未満の検証が
10,000の線形空間なのか
10,100の線形空間なのか
定義がずれているわけで、
AIやらニューロネットワーク的には、後者を選択するわけですよね？
52:卵の名無しさん
17/06/10 16:25:55.95 6oyq1OSX.net
明らかにn=10,000がn->無限と規定したいるから
おかしい話になるわけだ。
53:卵の名無しさん
17/06/16 15:01:52.10 B8W+f/XQ.net
お久しぶりです、スレ立てしてもメンテができず
メンテしてくれている人に感謝です。
政治ネタは下げたいのでご容赦を。
安倍政権のお友達の確率を二項分布でpとしたときに
明らかに3つ出れば有意な差が出る
でしょうね。
医学部を増やす、何の意味があるのでしょうか
竹中さん？
と思いますよ、今の医学系の経営リスクを増やすんでしょうな。
54:卵の名無しさん
17/06/16 15:07:56.72 B8W+f/XQ.net
>>50
さて先生のご不満はご承知の上で挑みますよ。
臨床統計における
独立したヒルベルト空間における独立したベクトルをどう設定するかでしょ？
主治医が男性か女性、あるいは患者さんにとって異性か同性なのか、
それを独立したベクトルとして採用するかどうか
が問題なんでしょうね。
臨床統計にとって、
ヒルベルト空間で何が独立した単位ベクトルなのか
もちろん、集めた母集団で主因子分析をするわけですが
「しゅいんし」とかな漢字変換に渡したら、「手淫し分析」
が一番に上がるぐらい自己満足感満載です。
55:卵の名無しさん
17/06/16 15:12:47.34 B8W+f/XQ.net
ここ数日の国会中継を見ていますが、
加世学園の問題は確率論的にベイズ更新を行って
安倍総理の
ここ数年の法案提出 → 国会での可決
確率pがが落ちているのを
覚悟すべきだと思います。
既にp自体が低下しているのに泥船に乗る確率を
落とすべきですよね。
ゲーム理論で表現される、
混合戦略に入りつつあるわけです。
56:卵の名無しさん
17/06/17 18:17:54.34 Ew/SZFmH.net
中元を配布したリストの提出を求められて税務署に提出。
税務署が「無作為」抽出（実は無作為抽出でなく作為抽出）して調査した５例中、中元を受け取ったのは０であったという。
それをもって税務署は中元は１００例全例虚偽であると認定した。
これはサンプリングに伴うバラつきだと主張して全例への課税を軽減したい。
どういう計算をすればいいか？
57:卵の名無しさん
17/06/18 16:18:56.05 nxMdZkja.net
# NN(100)個中当たりSS個、N(5)個サンプルしてS個当たりからSSを推測する。
foo <-function(SS,NN=100,N=5){
YY=c(rep(1,SS),rep(0,NN-SS))
return(sum(sample(YY,N)))
}
ss=0:100
S=0
k=10^4
SS=NULL
for(i in 1:k){
x=sapply(ss,foo)
SS=c(SS,which(x==S)-1)
}
hist(SS,freq=FALSE,col='lightblue', main='')
lines(density(SS),lwd=2,lty=3)
summary(SS)
MAP(SS)[1]
quantile(SS,c(0.025,0.50,0.975))
> summary(SS)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 4.00 10.00 13.56 20.00 89.00
> MAP(SS)[1]
x
1.547634
> quantile(SS,c(0.025,0.50,0.975))
2.5% 50% 97.5%
0 10 45
URLﾘﾝｸ(i.imgur.com)
58:卵の名無しさん
17/06/18 16:40:37.77 nxMdZkja.net
当たり確率を連続量のpとしてstanを走らせると似たような結果が得られた。
> print(fit, pars='p',probs=c(.025,.5,.975),digits=4)
Inference for Stan model: coin5.
4 chains, each with iter=100000; warmup=50000; thin=1;
post-warmup draws per chain=50000, total post-warmup draws=200000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
p 0.1429 0.0005 0.123 0.0042 0.1097 0.4569 68121 1
URLﾘﾝｸ(i.imgur.com)
59:卵の名無しさん
17/06/19 06:33:53.85 chw3mnRg.net
>>57
5C0*x^0*(1-x)^(5-0)の期待値を求めればいいので
N=5
S=0
integrate(function(x) choose(N,S)*x^S*(1-x)^(N-S),0,1)$value
[1] 0.1666667
シミュレーションと少し違うな。
60:卵の名無しさん
17/06/19 06:36:19.05 chw3mnRg.net
中元を送った期待値はNN/(N+1)でいいかな？
61:卵の名無しさん
17/06/22 15:04:26.25 StHWKWcJ.net
お久さです。はあ、
世田谷区でデング熱と診断された人が蚊に刺された
として、デング熱の発症確率は増えるんでしょうか？
根拠はあるのですか？ 行政的な出費をする根拠になるのですか？
62:卵の名無しさん
17/07/01 15:48:51.40 5TEeU/0W.net
お久しぶりですが、
AIの登場で臨床統計も変わってくると思いますよね。
AIなんてベイズ統計学そのもので検定主義の統計学とは異なるわけです。
死亡率と治療手段を独立とするベクトル空間における内積で議論するには
問題があるように思えますが。
個々の患者さんを説得するには明らかに分散が変わってくると思いますよね？
63:卵の名無しさん
17/07/01 20:17:24.70 UQXead44.net
>>62
ベイズ統計で進めていくときに
説明変数が独立かどうかの設定が難しいと思う。
どうやるのか知らないからではあるが。
模擬試験の判定って ロジスティック回帰の応用だろうけど
英語の点数と国語の点数って独立じゃないだろね。
最近は、ジェネラリストのための内科診断リファレンス
をつまみ読みして生半可な知識を強化している。
明日は雨でしょうから
明日の降水確率は**%
に昇華。
知識が正確になると却って混乱することも増えた。
バセドー病でも抗TPO抗体が陽性になることもあるとか。
64:卵の名無しさん
17/07/21 15:51:33.35 u9ETFd0W.net
ドクターG
URLﾘﾝｸ(www4.nhk.or.jp)
ありますよね。今でもNHK-Gで毎週放送しています。
新感覚！病名推理エンターテインメント番組「総合診療医ドクターG」。
病名を探り当てるまでの謎解きの面白さをスタジオで展開する！
あなたの症状も解き明かされるか！
というNHK的なクレジットがあるのですが、
患者さんの生死とかQOKがゲーム的に扱われている
のに疑問を感じます。
診断をどう扱うか、確率のヒルベルト空間では
ベイズ更新的
な発想で考えていますよね？
一つの事実で1,000-10,000ある診断ラベルのどれの
確率が変わるのか
という話です。
統計学・AI的思考の話。
研修医は診断AIに学ぶべき
だと思いませんか？
65:卵の名無しさん
17/07/21 16:01:59.58 u9ETFd0W.net
ある確率密度関数
P(x)
に関して
時間tの要素があるんじゃないか
P(x,t)と思いますよね。
2000年の評価と2017年の評価は違う
としたら、
2000年に治療を受けていた患者さんと
2017年に治療を受けていた患者さんは
別だ
と言いたいわけでしょ？
それで
サイコロを2回振るのか二つのサイコロを1回振るのか
の違いが分かると思いますよ。
66:卵の名無しさん
17/07/21 16:05:13.19 u9ETFd0W.net
今日、午後4時にNHK-Gにチャンネルを合わせると
白鵬の土俵入りが観られる
という偶然は偶然じゃなくて
それなりに
ポテンシャルを変更している
という話が量子論的・確率論的な話だよね
67:。
68:卵の名無しさん
17/07/21 16:08:29.61 u9ETFd0W.net
千秋楽に
2人の横綱が対決する
ことを事後確率にして
何人の横綱が必要か
という問題を解いているわけだね。
69:卵の名無しさん
17/08/17 20:34:13.43 XB8FFU5P.net
医者ならば、シリツ卒なら馬鹿である から
シリツ卒ならば、医者ならば馬鹿である が、導けるか？
という論理命題の問題を臨床適用ｗするとこうなる。
甲状腺癌ならば、未分化癌なら予後不良である から
未分化癌ならば、甲状腺癌なら予後不良である が、導けるか？
Rに真理表を使って解答させるス�
70:Nリプトは以下の通り。 # P->(Q->R) |- Q->(P->R) library('gtools') pm3=permutations(2,3,v=c(T,F),re=TRUE) colnames(pm3)=c('P','Q','R'); pm3 imply <- function(x,y) !(x&&!y) f <- function(P,Q,R) imply(imply(P,imply(Q,R)),imply(Q,imply(P,R))) f1 <- function(pm) f(pm[1],pm[2],pm[3]) result=logical(8) for(i in 1:8) result[i]=f1(pm3[i,]) cbind(pm3,result) > cbind(pm3,result) P Q R result [1,] FALSE FALSE FALSE TRUE [2,] FALSE FALSE TRUE TRUE [3,] FALSE TRUE FALSE TRUE [4,] FALSE TRUE TRUE TRUE [5,] TRUE FALSE FALSE TRUE [6,] TRUE FALSE TRUE TRUE [7,] TRUE TRUE FALSE TRUE [8,] TRUE TRUE TRUE TRUE
71:卵の名無しさん
17/09/15 02:08:09.73 qA2A3XQj.net
ここは臨床統計？
72:卵の名無しさん
17/09/15 02:58:26.93 zP5c6quI.net
>>69
>50できる？
73:卵の名無しさん
17/09/15 03:01:19.88 zP5c6quI.net
>>70
(3)若い女医の担当患者の成人率が、老いた男医の担当患者の生存率より高いと言えるか？
も追加。
74:卵の名無しさん
17/09/24 19:46:12.50 ZHgpf3LZ.net
「くじ引きが無作為である」という帰無仮説のもとで宝くじに当選する確率はとても低い（０．０５未満)。
宝くじに当選者がでたということはp<0.05のことが起こったので「くじ引きが無作為」という帰無仮説は棄却される。
正しい議論か？
75:卵の名無しさん
17/09/24 20:21:22.18 vQO8g+Lx.net
>>56
# NN(100)個中当たりSS個、N(5)個サンプルしてS個当たりからSSを推測する。
foo <-function(SS,NN=100,N=5){
YY=c(rep(1,SS),rep(0,NN-SS)) #YYは当たり(1)がSS個、はずれ(0)がNN-SS個
return(sum(sample(YY,N))) # YYからN個引いたときの当たりの個数を返す
}
sss=0:100 #くじ中のあたりの候補
S=0 #引いた当たりの数
k=10^4 # 繰り返し回数
SS=NULL #容れ子初期値
for(i in 1:k){
x=sapply(sss,foo) #あたりがくじ全体でsss個だった時N個引いたときのあたりの個数
SS=c(SS,which(x==S)-1) #引いたあたりの個数がS個(0個)だったときのくじ全体のあたりの個数で配列をつくる
}
hist(SS,freq=FALSE,col='lightblue', main='')
lines(density(SS))
summary(SS)
quantile(SS,c(0.025,0.50,0.975))
## 最頻値計算
MAP <- function(x) {
dens <- density(x)
mode_i <- which.max(dens$y) # densityの頂点となるindex
mode_x <- dens$x[mode_i]
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
MAP(SS)[1] #最頻値
76:卵の名無しさん
17/09/24 22:09:52.30 vQO8g+Lx.net
> summary(SS)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 4.00 10.00 13.56 20.00 85.00
> quantile(SS,c(0.025,0.50,0.975))
2.5% 50% 97.5%
0 10 45
> NN=100
> pp=0:NN #くじ全体のアタリの個数候補
> f <- function(x,NN=100,N=5) choose(NN-x,N)/choose(NN,N) #アタリがx個のときN個引いてすべてハズレの確率
> plot(pp,f(pp))
> sum(pp*f(pp)/sum(f(pp))) # f(pp)/sum(f(pp)):確率密度
[1] 13.57143
77:卵の名無しさん
17/09/24 22:53:55.42 vQO8g+Lx.net
# ド底辺特殊シリツ医大は1学年100人として全学生600人とする。このうち20人を調査したところ全員裏口入学であった。全学年での正規入学者数の期待値とその95%信頼区間を求めよ。
# 期待値＝26.36人 信頼区間0～101.06人
NN=600 ; N=20
pp=0:NN #くじ全体のアタリの個数候補
f <- function(x,NN=600,N=20) choose(NN-x,N)/choose(NN,N) #アタリがx個のときN個引いてすべてハズレの確率
plot(pp,f(pp))
sum(pp*f(pp)/sum(f(pp))) # f(pp)/sum(f(pp)):確率密度
binom::binom.confint(0,N)
(p=1-0.025^(1/N))
(1-p)^N ; p*NN
binom.test(0,N)$conf*NN
binom::binom.exact(0,N)
# NN個中アタリSS個、N個サンプルしてS個アタリからSSを推測する。
foo <-function(SS,NN=600,N=20){
YY=c(rep(1,SS),rep(0,NN-SS)) #YYはアタリ(1)がSS個、はずれ(0)がNN-SS個
return(sum(sample(YY,N))) # YYからN個引いたときのアタリの個数を返す
}
sss=0:NN #くじ中のアタリの候補
S=0 #引いたアタリの数
k=10^4 # 繰り返し回数
SS=NULL #容れ子初期値
for(i in 1:k){
x=sapply(sss,foo) #アタリがくじ全体でsss個だった時N個引いたときのアタリの個数
SS=c(SS,which(x==S)-1) #引いたアタリの個数がS個(0個)だったときのくじ全体のアタリの個数で配列をつくる
}
hist(SS,freq=FALSE,col='wheat', main='',xlab='正規合格者数') ; lines(density(SS))
summary(SS) ; quantile(SS,c(.025,.50,.975))
conf=binom::binom.confint(0,20)
cbind(conf['method'],round(conf['lower'],1),round(conf['upper']*NN,2))
78:卵の名無しさん
17/09/24 23:26:12.20 vQO8g+Lx.net
期待値 sum(pp*f(pp)/sum(f(pp)))
n=600
r=20
n
79:ｎ Σj*(n-j)Ｃr/nＣr ÷ Σ(n-i)Ｃr/nＣr j=0 i=0
80:卵の名無しさん
17/09/25 01:09:04.22 1DSpMtz1.net
## 全体.NN個中アタリSS個、.N個引いてr個アタリからSSを推測する。
atari <- function(.NN, .N, r, k=10^3){
f <-function(SS,NN=.NN,N=.N){
YY=c(rep(1,SS),rep(0,NN-SS)) #YYはアタリ(1)がSS個、はずれ(0)がNN-SS個
sum(sample(YY,N)) # YYからN個引いたときのアタリの個数を返す
}
sss=0:.NN # くじ中のアタリの候補
SS=NULL # 容れ子初期値
for(i in 1:k){ # k:繰り返し数
x=sapply(sss,f) # アタリがくじ全体でsss個だった時N個引いたときのアタリの個数
SS=c(SS,which(x==r)-1) #引いたアタリの個数がr個だったときのくじ全体のアタリの個数で配列をつくる
}
hist(SS,freq=FALSE,col='lightblue', main='', xlab='アタリ総数')
lines(density(SS))
print(summary(SS))
cat('\n')
print(quantile(SS,c(.025,.05,.50,.95,.975)))
invisible(SS)
}
atari(100,5,0)
atari(600,20,1)
### くじ全体のアタリの期待値
atari.e <- function(NN,N,r){ # 全体NN個のくじからN個引いてr個アタリ
pp=0:NN #くじ全体のアタリの個数候補
f <- function(x) choose(x,r)*choose(NN-x,N-r)/choose(NN,N) # 全体のアタリがx個のときN個引いてr個アタリの確率
plot(pp,f(pp))
sum(pp*f(pp)/sum(f(pp))) # Σpp：個数 * f(pp)/sum(f(pp)):確率密度
}
atari.e(100,5,0)
atari.e(600,20,1)
81:卵の名無しさん
17/09/25 06:39:09.92 1DSpMtz1.net
# N N
# Σ i*iＣr*(N-i)Ｃ(n-r)/NＣn ÷ ΣjＣr*(N-j)Ｃ(n-r)/NＣn
# i=0 j=0
# Σi*choose(i,r)*choose(N-i,n-r)/choose(N,n) ÷ Σchoose(j,r)*choose(N-j,n-r)/choose(N,n)
atari.ee <- function(N,n,r){
X=numeric() ; Y=numeric()
for(i in 1:(N+1)) X[i] = (i-1)*choose(i-1,r)*choose(N-(i-1),n-r)/choose(N,n)
for(j in 1:(N+1)) Y[j] = choose(j-1,r)*choose(N-(j-1),n-r)/choose(N,n)
sum(X)/sum(Y)
}
atari.ee(100,5,0)
atari.ee(600,20,1)
82:卵の名無しさん
17/09/25 06:40:55.63 1DSpMtz1.net
# N N
# Σ i*iＣr*(N-i)Ｃ(n-r)/NＣn ÷ ΣjＣr*(N-j)Ｃ(n-r)/NＣn
# i=0 j=0
# Σi*choose(i,r)*choose(N-i,n-r)/choose(N,n) ÷ Σchoose(j,r)*choose(N-j,n-r)/choose(N,n)
atari.ee <- function(N,n,r){
X=numeric() ; Y=numeric()
for(i in 1:(N+1)) X[i] = (i-1)*choose(i-1,r)*choose(N-(i-1),n-r)/choose(N,n)
for(j in 1:(N+1)) Y[j] = choose(j-1,r)*choose(N-(j-1),n-r)/choose(N,n)
sum(X)/sum(Y)
}
atari.ee(100,5,0)
atari.ee(600,20,1)
83:卵の名無しさん
17/09/25 07:58:57.17 1DSpMtz1.net
##
N=600
n=20
r=1
L <- function(x) choose(x,r)*choose(N-x,n-r)/choose(N,n) # 全体のアタリがx個のときn個引いてr個アタリの確率
L(0)
L(1)
xx=1:N
yy=sapply(xx,L)
which.max(yy)
plot(xx,yy)
plot(xx[25:35],yy[25:35])
L(29:31)
84:卵の名無しさん
17/09/25 09:02:50.99 CFjnZboG.net
N=600
n=20
r=1
L <- function(x) choose(x,r)*choose(N-x,n-r)/choose(N,n) # 全体のアタリがx個のときn個引いてr個アタリの確率
L(0) # = 0
L(1)
xx=1:N # yy[i]=L(xx[i])になるように1から始める
yy=sapply(xx,L)
which.max(yy) # 最頻値
plot(xx,yy)
plot(xx[25:35],yy[25:35])
L(29:31)
#chooseはΓ関数で内部計算なので整数でなくても計算する
curve(L(x),0,600)
optimise(L,c(0,600),maximum=TRUE)
85:卵の名無しさん
17/09/25 15:08:29.67 uD79ucO8.net
>>156
素数Aiを任意の数ｎ個集めて
n
Ｂ = Π Ai + 1
i=1
とすると
Ｂを割り切るAi以外の素数が存在することになり、ｎ個集めればn+1個以上の素数の存在が示せることになる。
例
2*3+1で素数7
3*5+1で素数2
3*7+１で素数2と11
86:卵の名無しさん
17/09/25 20:13:55.84 uD79ucO8.net
# N N
# Σ i*iＣr*(N-i)Ｃ(n-r)/NＣn ÷ ΣjＣr*(N-j)Ｃ(n-r)/NＣn
# i=0 j=0
# Σi*choose(i,r)*choose(N-i,n-r)/choose(N,n) ÷ Σchoose(j,r)*choose(N-j,n-r)/choose(N,n)
atari.m <- function(N,n,r){
X=numeric() ; Y=numeric()
for(i in 1:(N+1)) X[i] = (i-1)*choose(i-1,r)*choose(N-(i-1),n-r)/choose(N,n)
for(j in 1:(N+1)) Y[j] = choose(j-1,r)*choose(N-(j-1),n-r)/choose(N,n)
c(mean=sum(X)/sum(Y),mode=which.max(Y)-1)
}
atari.m(100,5,0)
atari.m(600,20,1)
87:卵の名無しさん
17/09/25 20:32:22.48 uD79ucO8.net
# N N
# Σ i*iＣr*(N-i)Ｃ(n-r)/NＣn ÷ ΣjＣr*(N-j)Ｃ(n-r)/NＣn
# i=0 j=0
# Σi*choose(i,r)*choose(N-i,n-r)/choose(N,n) ÷ Σchoose(j,r)*choose(N-j,n-r)/choose(N,n)
atari.mm <- function(N,n,r){
xx=0:N
c(mean = sum(xx*choose(xx,r)*choose(N-xx,n-r)/choose(N,n))/sum(choose(xx,r)*choose(N-xx,n-r)/choose(N,n)),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1)
}
atari.mm(100,5,0)
atari.mm(600,20,1)
88:卵の名無しさん
17/09/25 20:35:29.86 uD79ucO8.net
> atari.mm <- function(N,n,r){
+ xx=0:N
+ c(mean = sum(xx*choose(xx,r)*choose(N-xx,n-r)/choose(N,n))/sum(choose(xx,r)*choose(N-xx,n-r)/choose(N,n)),
+ mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1)
+ }
> atari.mm(100,5,0)
mean mode
13.57143 0.00000
> atari.mm(600,20,1)
mean mode
53.72727 30.00000
>
89:卵の名無しさん
17/09/25 20:38:01.19 uD79ucO8.net
# N N
# Σ i*iＣr*(N-i)Ｃ(n-r)/NＣn ÷ ΣjＣr*(N-j)Ｃ(n-r)/NＣn
# i=0 j=0
# Σi*choose(i,r)*choose(N-i,n-r)/choose(N,n) ÷ Σchoose(j,r)*choose(N-j,n-r)/choose(N,n)
atari.mm <- function(N,n,r) c(mean = sum(0:N*choose(0:N,r)*choose(N-0:N,n-r)/choose(N,n))/sum(choose(0:N,r)*choose(N-0:N,n-r)/choose(N,n)),mode = which.max(choose(0:N,r)*choose(N-0:N,n-r))-1)
atari.mm(100,5,0)
atari.mm(600,20,1)
90:卵の名無しさん
17/09/25 20:48:49.89 uD79ucO8.net
>83が可読性が一番いいな
91:卵の名無しさん
17/09/26 10:38:08.14 dTnYYVrv.net
# 全体.N個中アタリS個、1個ずつクジを引いて当たったらやめる. r個めがアタリからSを推測する。
lottery <- function(.N,.r,k=10^3){
f <-function(S,N=.N){
y=c(rep(1,S),rep(0,N-S)) # yはアタリ(1)がS個、はずれ(0)がN-S個
Y=sample(y,N) # 並べ替えた配列を返す
return(Y)
}
# 初めて0以外の数が現れた位置を返す
g <- function(y){
n=length(y)
for(i in 1:n){
if(!y[i]) i=i+1
else break
}
return(i)
}
xx=0:.N # くじ中のアタリ数の候補
SS=NULL # 容れ子初期値
for(i in 1:k){
M=t(sapply(xx,f)) # アタリ0~.N個の並べ替え済サンプル行列M
# r個めが初めて当たりだったときの個数（=行番-1)
for(j in 1:.N){
if(g(M[j,])==.r) SS=c(SS,j-1) #r個めで初めてのアタリだったときのくじ全体のアタリの個数で配列をつくる
}
}
hist(SS,freq=FALSE,main='', xlab='Wins',col='skyblue')
# lines(density(SS))
print(summary(SS))
cat('\n')
print(quantile(SS,c(.025,.05,.50,.95,.975)))
invisible(SS)
}
92:卵の名無しさん
17/09/26 10:38:34.73 dTnYYVrv.net
##
lottery.mm <- function(.N,.r){
pmf <- function(S,N=.N,r=.r){ # .N個中S個のアタリのとき.r個めで初めてあたる確率は
choose(N-r,S-1)/choose(N,S) #.r+1以後の配列/すべての配列
}
ss=0:.N
pdf=pmf(ss)/sum(pmf(ss)) # 確率密度関数に変換
plot(ss,pdf)
c(mean=sum(ss*pdf),mode=which.max(pdf)-1)
}
lottery.mm(100,50)
93:卵の名無しさん
17/09/26 15:09:39.24 dTnYYVrv.net
lottery.mci <- function(.N,.r){ # 期待値と信頼区間を求める
pmf <- function(S,N=.N,r=.r){ # .N個中S個のアタリのとき.r個めで初めてあたる確率は
choose(N-r,S-1)/choose(N,S) #.r+1以後の配列/すべての配列
}
ss=0:.N
pdf=pmf(ss)/sum(pmf(ss)) # 確率密度関数に変換
E=sum(ss*pdf)
SSS=sample(0:.N,10^5*.N,replace=TRUE,prob=pdf) # 重みをつけて繰り返しサンプリング
UL=quantile(SSS,c(.025,.975))
c(E,UL)
}
lottery.mci(100,5)
94:卵の名無しさん
17/09/26 19:18:43.02 VKfC0TE5.net
> lottery.mci <- function(.N,.r){ #
+ pmf <- function(S,N=.N,r=.r){ #
+ choose(N-r,S-1)/choose(N,S) #
+ }
+ ss=0:.N
+ pdf=pmf(ss)/sum(pmf(ss)) #
+ E=sum(ss*pdf)
+ SSS=sample(0:.N,10^5*.N,replace=TRUE,prob=pdf) #
+ UL=quantile(SSS,c(.025,.975))
+ c(mean=E,UL)
+ }
> lottery.mci(100,1)
mean 2.5% 97.5%
67 16 99
95:卵の名無しさん
17/09/26 20:12:10.07 VKfC0TE5.net
lottery.mci <- function(.N,.r){
pmf <- function(S,N=.N,r=.r){
choose(N-r,S-1)/choose(N,S)
}
ss=0:.N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:.N,10^5*.N,replace=TRUE,prob=pdf)
UL=quantile(SSS,c(.025,.975))
c(E,UL)
}
.N=100
rr=1:.N
MLU=sapply(rr, function(r) lottery.mci(.N,r))
round(t(MLU))
96:卵の名無しさん
17/09/26 20:25:43.72 VKfC0TE5.net
> round(t(MLU))
2.5% 97.5%
[1,] 67 16 99
[2,] 50 9 91
[3,] 40 7 80
[4,] 33 5 71
[5,] 28 4 63
[6,] 25 4 57
[7,] 22 3 51
[8,] 19 3 47
[9,] 18 2 43
[10,] 16 2 40
[11,] 15 2 37
[12,] 14 2 34
[13,] 13 2 32
[14,] 12 2 30
[15,] 11 2 28
[16,] 10 1 27
[17,] 10 1 25
[18,] 9 1 24
[19,] 9 1 23
[20,] 8 1 22
[21,] 8 1 21
[22,] 8 1 20
[23,] 7 1 19
[24,] 7 1 18
[25,] 7 1 17
[26,] 6 1 17
[27,] 6 1 16
[28,] 6 1 15
[29,] 6 1 15
[30,] 5 1 14
97:卵の名無しさん
17/09/27 08:43:44.49 lL+6eziV.net
atari.Q <- function(N,n,r,k=10^3){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
print(quantile(SS,c(0.025,0.5,0.975)))
c(mean = sum(xx*pdf),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1)
}
atari.Q(100,10,2)
98:卵の名無しさん
17/09/27 08:44:44.68 lL+6eziV.net
> atari.Q(100,10,2)
2.5% 50% 97.5%
6 23 50
mean mode
24.5 20.0
>
99:卵の名無しさん
17/09/27 10:05:45.43 UOZ+ZQoW.net
> lottery.mm <- function(.N,.r){
+ pmf <- function(S,N=.N,r=.r){
+ choose(N-r,S-1)/choose(N,S)
+ }
+ ss=0:.N
+ pdf=pmf(ss)/sum(pmf(ss))
+ c(mean=sum(ss*pdf),mode=which.max(pdf)-1)
+ }
> lottery.mm(100,5)
mean mode
28.14286 20.00000
100:卵の名無しさん
17/09/27 16:26:42.51 By0OihfH.net
N=100 ; n=5 ; r=0 ; m=10
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
print(c(mean = sum(xx*pdf),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
cat('\n')
print(quantile(SS,c(0.025,0.5,0.975)))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
print(c(mean=E,mod=which.max(pdf)))
cat('\n')
print(quantile(SSS,c(.025,.975)))
invisible(SSS)
}
atari.Q(N,n,r)
lottery.mci(N,m)
101:卵の名無しさん
17/09/27 16:27:31.48 By0OihfH.net
>>97
シリツ医大進学予備校"どていへん予備校"では５年連続不合格、
別のシリツ医大進学予備校"うらぐち予備校"では１０年めで初めて合格者がでたという実績があるとき、
どていへん予備校 と うらぐち予備校ではどちらが合格可能性が高いと言えるか？
102:innuendo
17/09/27 16:37:18.00 By0OihfH.net
巨乳女子大100人と桃尻女子大100人のどちらが誘いに応じやすいか検討してみる。
巨乳女子大では10人に声をかけて2人が誘いに応じた。
桃尻女子大では無作為順に声をかけていったら5人めが誘いに応じた。
どちらが誘いに応じやすいといえるか？
URLﾘﾝｸ(i.imgur.com)
> lottery.mci(100,5)
mean mod
28.14286 21.00000
2.5% 97.5%
4 63
> atari.Q(100,10,2)
mean mode
24.5 20.0
2.5% 50% 97.5%
6 23 50
>
103:卵の名無しさん
17/09/27 17:32:04.45 By0OihfH.net
N=100 ; n=5 ; r=0 ; m=10
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
print(c(mean = sum(xx*pdf),
mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
cat('\n')
print(quantile(SS,c(0.025,0.5,0.975)))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
print(c(mean=E,mod=which.max(pdf)-1))
cat('\n')
print(quantile(SSS,c(.025,.975)))
invisible(SSS)
104:卵の名無しさん
17/09/27 17:40:20.17 By0OihfH.net
N=100 ; n=10 ; r=2 ; m=5 # バグ修正版
atari.Q <- function(N,n,r,k=10^5){
xx=0:N
pmf=choose(xx,r)*choose(N-xx,n-r)/choose(N,n)
pdf=pmf/sum(pmf)
SS=sample(xx,N*k,replace=TRUE,prob=pdf)
# print(quantile(SS,c(0.025,0.5,0.975)))
# cat('\n')
# print(c(mean = sum(xx*pdf),
# mode = which.max(choose(xx,r)*choose(N-xx,n-r))-1))
invisible(SS)
}
lottery.mci <- function(N,m,k=10^5){
pmf <- function(S,.N=N,.r=m){
choose(.N-.r,S-1)/choose(.N,S)
}
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
E=sum(ss*pdf)
SSS=sample(0:N,N*k,replace=TRUE,prob=pdf)
UL=quantile(SSS,c(.025,.975))
# print(c(E,UL))
invisible(SSS)
}
Brown=sum(choose(0:N,r)*choose(N-0:N,n-r)/choose(N,n))
curve(choose(x,r)*choose(N-x,n-r)/choose(N,n)/Brown,0,100,col='brown',xlab='尻軽数',ylab='')
Pink=sum(choose(N-m,0:N-1)/choose(N,0:N))
curve(choose(N-m,x-1)/choose(N,x)/Pink,0,100, add=TRUE,col='pink',lwd=2)
legend('topright',bty='n',legend=c('巨乳女子','桃尻女子'),lwd=1:2,col=c('brown','pink'))
105:innuendo
17/09/27 17:59:40.63 By0OihfH.net
URLﾘﾝｸ(i.imgur.com)
平均値と信頼区間のどちらを選択するかだな。
106:卵の名無しさん
17/09/27 20:59:53.74 By0OihfH.net
options(scipen = 100)
lottery.s <- function(m,N=100){
pmf <- function(S,.N=N,.m=m){
choose(.N-.m,S-1)/choose(.N,S)
}
bb=0:N # numbers of possible bitches
pdf=pmf(bb)/sum(pmf(bb))
SS=sample(0:N,replace=TRUE,prob=pdf) # sample as the same length
return(SS)
}
lottery.s(8)
foo <- function(y,al='l'){ # one-sided mo > ky
f <- function(x){
mo=lottery.s(5)
ky=lottery.s(x)
t.test(ky,mo,alt=al)$p.value
}
welch=replicate(10^2,f(y))
mean(welch)
}
yy=5:10
pp=sapply(yy,foo)
plot(yy,pp)
abline(h=0.05,lty=3)
107:卵の名無しさん
17/09/27 21:27:00.63 By0OihfH.net
foo <- function(y,fn=t.test,al='l'){ # one-sided mo > ky
f <- function(x){
mo=lottery.s(5)
ky=lottery.s(x)
fn(ky,mo,alt=al)$p.value
}
p.values=replicate(10^2,f(y))
mean(p.values)
}
yy=5:10
# どの検定でも結果は変わらず
pp=sapply(yy,function(x) foo(x,fn=t.test))
pp=sapply(yy,function(x) foo(x,fn=lawstat::brunner.munzel.test))
pp=sapply(yy,function(x) foo(x,fn=wilcox.test))
PP=sapply(yy,function(x) foo(x,fn=perm::permTS))
plot(yy,pp)
abline(h=0.05,lty=3)
108:卵の名無しさん
17/09/27 23:32:37.60 By0OihfH.net
> lottery.e <- function(N,m){
+ pmf <- function(S,.N=N,.m=m){
+ choose(.N-.m,S-1)/choose(.N,S)
+ }
+ ss=0:N
+ pdf=pmf(ss)/sum(pmf(ss))
+ sum(ss*pdf)
+ }
>
> lottery.e(100,10)
[1] 16
> lottery.e(50,5)
[1] 13.85714
>
> A=100;a=10
> B=50;b=5
>
> prop.test(c(lottery.e(A,a),lottery.e(B,b)),c(A,B))
2-sample test for equality of proportions with continuity correction
data: c(lottery.e(A, a), lottery.e(B, b)) out of c(A, B)
X-squared = 2.1814, df = 1, p-value = 0.1397
alternative hypothesis: two.sided
95 percent confidence interval:
-0.27551116 0.04122545
sample estimates:
prop 1 prop 2
0.1600000 0.2771429
109:卵の名無しさん
17/09/28 00:10:25.54 bltMOtIo.net
sillygal <- function(a,b,A,B){
lottery.e <- function(N,m){
pmf <- function(S,.N=N,.m=m) choose(.N-.m,S-1)/choose(.N,S)
ss=0:N
pdf=pmf(ss)/sum(pmf(ss))
sum(ss*pdf)
}
prop.test(c(lottery.e(A,a),lottery.e(B,b)),c(A,B))
}
sillygal(5,10,50,100)
sillygal (5,7,100,100)
110:卵の名無しさん
17/09/29 12:14:08.07 aSJfCsJ7.net
1 st 2 nd 3 rd 4 th と表示させるための1行スクリプト
th=ifelse(0<r&&r<4,c('st','nd','rd')[which(r==1:3)],'th')
111:卵の名無しさん
17/09/29 19:11:58.46 JC0gA3LG.net
N=100
plot(0:N/N,dbinom(0:N,N,p=(0:N)/N),pch='+',
xlab='治癒確率',ylab='治癒確率通り治癒する確率')
yy=dbinom(0:N,N,p=(0:N)/N)
summary(yy[2:100])
hist(yy[2:100],col='tomato')
112:卵の名無しさん
17/09/30 18:41:04.08 ZzasON2B.net
# 元利均等返済 毎月の返済額は一定
.D0=10^8 # 借入額
.ra=0.01 # 年利
.r=.ra/12 # 月利
.n=120 # 返済月数
L <- function(x,D0=.D0,r=.r,n=.n){ # x:毎月の返済額
D=numeric(n)
D[1]=D0*(1+r) - x
for(i in 1:(n-1)) D[i+1] <- D[i]*(1+r) -x
return(D[n]) # n回返済後の借入金残高を返す
}
# L(.n)=0になる毎月の返済額 x を求める。
(m=uniroot(L,c(.D0/.n,.D0*(1+.n*.r)))$root)
m*.n # 総返済額
##### 元利均等返済 毎月の返済額と総返済額
loan <- function(.D0,.ra,.n){
A0=.D0
r=.ra/12
N=.n
m=A0*r*(1+r)^N/( (1+r)^N- 1)
c(Monthly=floor(m),Total=floor(m)*N)
}
loan(23700000,0.03,60)
loan(10^8,0.01,120)
113:卵の名無しさん
17/09/30 18:41:31.88 ZzasON2B.net
## 元利均等返済 nケ月めの残高
Debt <- function(n,.D0,.ra,.n){
A0=.D0
r=.ra/12
N=.n
s=1+r
A0*(1- (1-s^n)/(1-s^N))
}
Debt(1:60,23700000,0.03,60)
##
Loan <- function(.D0,.ra,.n){
L <- function(x,D0=.D0,r=.ra/12,n=.n){
D=numeric(n)
D[1]=D0*(1+r) - x
for(i in 1:(n-1)) D[i+1] <- D[i]*(1+r) -x
return(D[n])
}
m=uniroot(L,c(.D0/.n,.D0*(1+.n*.ra/12)))$root
c(Monthly=floor(m),Total=floor(m)*.n)
}
114:卵の名無しさん
17/10/01 12:38:58.98 NynEuCM1.net
foo <- function(N){
#N=7 # 元数
I=1:N
.Y=sample(I,N) ; .Y # 1,2,..,Nを.Y[1],Y.[2],..,.Y[N]に置換
tikan <- function(X,Y=.Y){ # X:置換される配列
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z) # 置換後の配列を返す
}
## 循環したらやめる
z=matrix(1:N,nrow=1)
z[1,]=tikan(I)
for(i in 1:N^2){
if(!prod(z[i,]==I)){
z=rbind(z,tikan(z[i,]))
}
else break
}
print(z)
print(rbind(I,z[1,]))
return(nrow(z))
}
115:卵の名無しさん
17/10/01 12:50:45.81 NynEuCM1.net
junkan <- function(.Y){ # 何回目で巡回したか返す、巡回しなければ０
N=length(.Y) # 元数
I=1:N
#.Y=sample(I,N) ; .Y # 1,2,..,Nを.Y[1],Y.[2],..,.Y[N]に置換
tikan <- function(X,Y=.Y){ # X:置換される配列
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z) # 置換後の配列を返す
}
## 循環したらやめる
z=matrix(1:N,nrow=1)
z[1,]=tikan(I)
for(i in 1:N^2){
if(!prod(z[i,]==I)){
z=rbind(z,tikan(z[i,]))
}
else break
}
# print(rbind(I,.Y))
# print(z)
return(ifelse(i!=N^2,i,ifelse(prod(z[i,]==I),i,0)))
}
116:卵の名無しさん
17/10/01 12:51:32.13 NynEuCM1.net
#
library(gtools)
N=7
PM=permutations(N,N,v=1:N)
jn=numeric(N)
n=nrow(PM)
for(j in 1:n) jn[j] <- junkan(PM[j,])
table(jn)
which(jn==0)
> table(jn)
jn
1 2 3 4 5 6 7 10 12
1 231 350 840 504 1470 720 504 420
> which(jn==0)
integer(0)
117:卵の名無しさん
17/10/01 17:13:29.79 NynEuCM1.net
test,mt1,st1,mt2,st2,mc1,sc1,mc2,sc2
開眼片足立ち左秒,16.02,7.73,18.99,8.12,18.92,9.84,18.83,9.92
開眼片足立ち右秒,18.48,8.29,21.23,8.14,19.88,9.36,19.75,9.35
左握力kg,18.72,4.22,19.67,4.03,19.83,4.63,19.73,4.68
右握力kg,20.89,4.13,21.7,3.86,21.21,4.37,20.92,4.28
FRcm,22.36,6.08,24.83,5.89,23.51,4.97,22.45,4.67
10m歩行速度秒,6.61,0.75,5.95,0.67,6.53,1.18,6.65,1.32
立位体前屈cm,6.90,10.64,12.33,6.41,9.25,4.90,9.00,5.48
左片足立ち振り回,8.40,4.76,9.93,4.24,8.96,6.28,8.36,5.85
右片足立ち振り回,9.30,5.20,11.16,5.34,9.93,7.38,9.80,7.27
nt=30
nc=30
(d=read.csv('clipboard'))
mt2=d[,'mt2']
st2=d[,'st2']
mc2=d[,'mc2']
sc2=d[,'sc2']
# t検定(生データなし,等分散不問)
Welch.test=function(n1,n2,m1,m2,sd1,sd2){
T=(m1-m2)/sqrt(sd1^2/n1+sd2^2/n2)
df=(sd1^2/n1+sd2^2/n2)^2 / (sd1^4/n1^2/(n1-1)+sd2^4/n2^2/(n2-1))
p.value=2*pt(abs(T),df,lower.tail = FALSE)
return(p.value)
}
cbind(d['test'],Welch.test(nt,nc,mt2,mc2,st2,sc2))
118:卵の名無しさん
17/10/03 10:21:03.41 D0GaF7Tg.net
# 百発百中
Dbinom <- function(s,n,p){
choose(n,s)*p^s*(1-p)^(n-s)
}
Dbinom(100,100,0.95)
Dbinom(100,100,0.99)
s=99;n=100 # 百発99中
curve(Dbinom(s,n,x),ylab='',xlab='p')
AUC=integrate(function(p) Dbinom(s,n,p),0,1)$value
integrate(function(x) x*Dbinom(s,n,x)/AUC,0,1)$value
s=100;n=100 # 百発百中
curve(Dbinom(s,n,x),ylab='',xlab='p',0.9,1)
AUC=integrate(function(p) Dbinom(s,n,p),0,1)$value
integrate(function(x) x*Dbinom(s,n,x)/AUC,0,1)$value # Expected value
curve(Dbinom(s,n,x)/AUC)
curve(Dbinom(s,n,x)/AUC,0.9,1)
integrate(function(x) Dbinom(s,n,x)/AUC, 0,1)
pdf <- function(x) Dbinom(s,n,x)/AUC
p99 <- uniroot(function(x,u0=0.99) integrate(pdf,x,1)$value - u0,c(0,1))$root ; p99
N=10^7
hits=rbinom(N,n,p99)
hist(hits)
mean(hits==100)
119:卵の名無しさん
17/10/03 14:56:19.98 D0GaF7Tg.net
# ｎ発ｎ中 CL=0.99
# PMF(x)=nCn*x^n*(1-x)^(n-n)=x^n
# AUC=integrate(PMF,0,1)=[x^(n+1)/(n+1)]1,0 = 1/(n+1)
# PDF=PMF/AUC=(n+1)x^n
# CDF(u)=integrate(pdf(x),0,u)=[x^(n+1)]0,u = u^(n+1)
# (1-CL)^(1/(n+1)) # lower.limit
# exp(log(1-CL)/(n+1))
#
# Ex=integrate(x*PDF,0,1) = integrate((n+1)*x^(n+1),0,1)
# = (n+1)/(n+2)
#
sniper <- function(n,CL=0.99){
c(Lower=(1-CL)^(1/(n+1)), Expected=(n+1)/(n+2), Upper=1)
}
sniper(100,0.99)
nn=1:100
plot(nn,sapply(nn,function(x) sniper(x,CL=0.99)[1]),
xlab='ｎ発ｎ中',ylab='99%信頼区間下限値',pch=19)
120:卵の名無しさん
17/10/03 14:56:47.23 D0GaF7Tg.net
##
N=10^3
n=100
pn=sniper(n,0.99)[1] ; pn
hits=rbinom(N,n,pn)
hist(hits,freq=FALSE, main='')
mean(hits==n)
k=10^5
M=numeric(k)
foo <-function(){
hits=rbinom(N,n,pn)
mean(hits==n)
}
for(i in 1:k) M[i] <- foo()
summary(M)
hist(M,freq=FALSE,main='',xlab='全発命中割合')
121:卵の名無しさん
17/10/04 11:34:22.68 Fg8f7jBH.net
2-way Contingency Table Analysis URLﾘﾝｸ(statpages.info)
Epi::twoby2
122:卵の名無しさん
17/10/04 13:13:51.18 Fg8f7jBH.net
# logRR±1.96*sqrt(1/a-1/(a+b)+1/c-1/(c+d))
# logOR±1.96*sqrt(1/a+1/b+1/c+1/d)
HRCI <- function(a,b,c,d,CL=0.95){
Z=qnorm(1-(1-CL)/2)
HR=(a/(a+b)) /(c/(c+d))
se=sqrt(1/a-1/(a+b)+1/c-1/(c+d))
lwr=HR*exp(-Z*se)
upr=HR*exp( Z*se)
c(lwr,HR,upr)
}
123:innuendo
17/10/04 20:50:12.19 Fg8f7jBH.net
To calculate M1, the relative risks in each of the six studies were combined using a random effects meta-analysis to give a point estimate of 0.361 for the relative risk with a confidence interval of (0.248, 0.527).
The 95% confidence interval upper bound of 0.527 represents a 47% risk reduction, which translates into a risk increase of about 90% from not being on warfarin (1/0.527 = 1.898)
(i.e., what would be seen if the test drug had no effect). Thus, M1 (in terms of the hazard ratio favoring the control to be ruled out) is 1.898.
124:innuendo
17/10/04 21:03:03.78 Fg8f7jBH.net
The clinical margin M2 representing the largest acceptable level of inferiority was therefore set at 50% of M1.
M2 was calculated on the log hazard scale as 1.378, so that NI would be demonstrated provided the upper bound for the 95% confidence interval for C vs. T < 1.378.
exp(0.5*log(1/0.527))=1.378
125:innuendo
17/10/05 11:09:02.87 56lc+jur.net
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
126:卵の名無しさん
17/10/05 13:06:15.38 8CwGr/wg.net
The historical trials assessed Ch ? P, the effect of the active control compared to placebo.
The NI trial assesses T ? Cn , the effect of the test drug compared to active control. If the outcome for the active control is constant across studies,
then Cn = Ch, and the sum T ? Cn + Ch ? P = T ? P represents the effect of the test drug compared to placebo.
127:卵の名無しさん
17/10/05 15:45:07.34 8CwGr/wg.net
# Non-Inferiority Trial
# Fixed Margin Approach
M1M2 <- function(ub){ # ub :upper border of c.i. by meta-analysis
c(M1=1/ub,M2=exp(-log(ub)/2))
}
M1M2(.527)
## Synthesis Method risk ratio
NIs <- function(rr21,se21,rr10,se10){ # 0:placebo, 1:active control 2:test drug
se=sqrt(se21^2+(se10/2)^2) # log_se=sqrt(log(1/a-1/(a+b)+1/c-1/c+d)))
z=(log(rr21)+log(rr10)/2)/se
p=pnorm(-abs(z))
c(p.value=p,z=z)
}
NIs(1.39,0.216,0.361,0.154)
## fixed margin test
NIf <-function(rr21,se21,rr10,se10){
se=se21+se10/2
z=(log(rr21)+log(rr10)/2)/se
p=pnorm(-abs(z))
c(p.value=p,z=z)
}
NIf(1.39,0.216,0.361,0.154)
128:卵の名無しさん
17/10/07 05:54:48.32 427btajj.net
G15 <- function(.N, .n, r, k=10^3,...){
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=0:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1)
}
summary(SS)
}
G(10000,100,100)
129:卵の名無しさん
17/10/07 06:05:26.39 427btajj.net
G13 <- function(N,n,r){
pp=0:N
f <- function(x) choose(x,r)*choose(N-x,n-r)/choose(N,n)
sum(pp*f(pp)/sum(f(pp)))
}
G13(10000,100,100)
G13(10000,10,10)
G13(10000,1,1)
130:卵の名無しさん
17/10/07 16:41:37.63 jwh+lzlc.net
Go13 <- function(.N, .n, r, k=10^3){
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
summary(SS)
}
Go13(100,1,1)
Go13(1000,10,10)
Go13(10000,100,100)
131:卵の名無しさん
17/10/07 19:50:28.87 427btajj.net
Na=1000 ; Nb=1000
pa=0.1 ; pb=0.01
Wa=100 ; Wb=1000
A=c(rep(1,Na*pa),rep(0,Na-Na*pa)) ; A=sample(A)
B=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb)) ; B=sample(B)
n=100
Get <- function(a){
sum(sample(A,a))*Wa + sum(sample(B,n-a))*Wb
}
aa=0:n
k=10^4
MAX=replicate(k,which.max(sapply(aa,Get)))
summary(MAX)
hist(MAX,freq=FALSE)
lines(density(MAX),lwd=2)
132:卵の名無しさん
17/10/08 06:36:02.71 HZvOcnfD.net
>>128
パチンコのスレにこういう記載があった。
>バカは何故かパチンコで勝てると思い込んでいる
>バカは本来なら勝てるのに遠隔や不正で負けていると思い込んでいる
>バカは10分の1で10円が当たるクジより1万分の1で8000円が当たるクジの方が儲かると思ってる
>算数すらできないバカが必死に守って支えてきたのがパチンコ業界
これを読んでこんな問題を考えてみた。
宝くじAは1000本中100本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
どちらも売り出し価格は同じなので100本買うことにする。
どちらを何本買うときに賞金の期待値が最大になるか、シミュレーションしてみよ。
133:卵の名無しさん
17/10/08 12:43:27.27 jiwdTQ4W.net
MAP <- function(x) { # モード値を計算
dens <- density(x)
mode_i <- which.max(dens$y)
mode_x <- dens$x[mode_i]
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
Na=1000 ; Nb=1000
pa=0.1 ; pb=0.01
Wa=100 ; Wb=1000
n=100
A0=c(rep(1,Na*pa),rep(0,Na-Na*pa))
B0=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb))
Get <- function(a){ #Aをa本買ったときの賞金合計
sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb
}
takarakuji <- function(k=10^3,Histogram=FALSE,Print=TRUE,...){
MAX=replicate(k,which.max(sapply(aa,Get))-1)
if(Histogram) {hist(MAX,freq=FALSE,...) ; lines(density(MAX),lwd=2)}
if(Print)print(round(c(mean=mean(MAX),sd=sd(MAX),mode=MAP(MAX)[1]),2))
invisible(c(mean=mean(MAX),sd=sd(MAX),mode=MAP(MAX)[1]))
}
takarakuji(10^5,Histogram = TRUE,col='lightblue')
134:卵の名無しさん
17/10/08 18:38:07.36 jiwdTQ4W.net
Takarakuji <- function(k=10^2,.Na=1000,.Nb=1000,.pa=0.1,.pb=0.01,
.n=100,.Wa=100,.Wb=1000,Histogram=TRUE,Print=TRUE,...){
Na=.Na ; Nb=.Nb
pa=.pa ; pb=.pb
Wa=.Wa ; Wb=.Wb
n=.n
aa=0:n
A0=c(rep(1,Na*pa),rep(0,Na-Na*pa))
B0=c(rep(1,Nb*pb),rep(0,Nb-Nb*pb))
Get <- function(a)sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb
MAX=replicate(k,which.max(sapply(aa,Get))-1)
tbl_MAX=table(MAX)
MODE=tbl_MAX[which.max(tbl_MAX)]
mode=as.numeric(names(MODE))
if(Histogram) {
hist(MAX,freq=FALSE,main=paste('Na=',Na,'Nb=',Nb),...)
lines(density(MAX),lwd=1)
}
if(Print){
print(round(c(mean=mean(MAX),sd=sd(MAX),mode=mode),2))
}
invisible(c(mean=mean(MAX),sd=sd(MAX),mode=mode))
}
dev.off()
par(mfrow=c(2,2))
Takarakuji(.Na=10^3,.Nb=10^3,col='lightblue')
Takarakuji(.Na=10^5,.Nb=10^5,col='wheat')
Takarakuji(.Na=10^7,.Nb=10^7,col='maroon')
Takarakuji(.Na=10^8,.Nb=10^8,col='gray')
135:卵の名無しさん
17/10/09 08:02:11.41 Wked2yci.net
# 10分の1で10万円が当たるクジAと1万分の1で8000万円が当たるクジの方Bある
# クジ1本の値段は前者は2万円,後者は4万円とする.
# 購入金額100万円としてどのように買うとき賞金の期待値が最大になるか？
pa=0.1 #Aの当選率
Wa=10 #Aの賞金
ca=2 #Aのコスト
pb=1/10000
Wb=8000
cb=4
Cash=100 #購入額
# a,bは各々A,Bの購入数 2*a+4*b=100 a=(Cash-4*b)/2
award <- function(b){
a=(100-4*b)/2
sum(rbinom(a, 1, pa) * Wa) + sum(rbinom(b, 1, pb) * Wb)
}
bb=0:floor(Cash/cb)
AWD <- function(){
Award=sapply(bb,award)
c(which.max(Award)-1,max(Award))
}
k=10^6 ; res=replicate(k,AWD())
Max.B=res[1,]
hist(Max.B,freq=FALSE,breaks=20,col='lightblue')
lines(density(Max.B),lwd=1)
summary(Max.B)
Mode(Max.B)
Max.AWD=res[2,]
hist(Max.AWD,col='wheat',breaks=100)
summary(Max.AWD)
Mode(Max.AWD)
136:卵の名無しさん
17/10/09 11:07:15.42 Wked2yci.net<
137:> Rstudioのdplyrとtidyr の Cheetsheet https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf のデータは EDAWR というパッケージにあるらしいのだが、R3.4.2だと Warning in install.packages : package ‘EDAWR’ is not available (for R version 3.4.2) と警告がでた。 install.packages('devtools') devtools::install_github("rstudio/EDAWR") で警告を回避してインストールできた。 > cases country 2011 2012 2013 1 FR 7000 6900 7000 2 DE 5800 6000 6200 3 US 15000 14000 13000 > tidyr::gather(cases,'year','n',2:4) country year n 1 FR 2011 7000 2 DE 2011 5800 3 US 2011 15000 4 FR 2012 6900 5 DE 2012 6000 6 US 2012 14000 7 FR 2013 7000 8 DE 2013 6200 9 US 2013 13000 と動作してくれた。
138:卵の名無しさん
17/10/09 14:34:35.08 Wked2yci.net
library(tidyr)
library(EDAWR)
cases
cases %>% gather('year','n',2:4)
cases %>% gather('year','n',2:4) %>% spread(year,n) %>% arrange(country)
cases %>% arrange(country)
pollution
pollution %>% spread(size, amount)
pollution %>% spread(size, amount) %>% gather('size','amount',2:3) %>% arrange(desc(city))
pollution
library(reshape2)
cases
cases %>% melt(value.name ='cases')
pollution
pollution %>% acast(city ~ size)
139:innuendo
17/10/09 15:13:49.58 N8CcYMr4.net
>>127
URLﾘﾝｸ(imagizer.imageshack.com)
140:卵の名無しさん
17/10/09 15:39:38.20 Wked2yci.net
>>135
URLﾘﾝｸ(i.imgur.com)
141:卵の名無しさん
17/10/09 15:40:36.46 Wked2yci.net
> summary(g13)
Min. 1st Qu. Median Mean 3rd Qu. Max.
8992 9866 9933 9903 9972 10000
> which.max(table(g13))
10000
714
142:卵の名無しさん
17/10/09 19:26:11.72 Wked2yci.net
pa=0.1 ; pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
k=10^3 # シミュレーション回数
# Aをi買ったときの賞金総額
a <- function(i){sum(rbinom(i, 1, pa) * Wa) + sum(rbinom(100-i, 1, pb) * Wb)}
which.max2 <- function(x){which(x==max(x))-1}
# 最大賞金総額になるAの購入数（複数可）の配列を返す
b <- function(){
b1=NULL
b1=c(b1,which.max2
143:(sapply(0:n, a))) return(b1) } b2=NULL for(j in 1:k){ b2=c(b2,b()) } hist(b2,freq=FALSE,breaks=20,col=sample(colors(),2),ann=FALSE) lines(density(b2),lwd=1) summary(b2) which.max(table(b2))
144:卵の名無しさん
17/10/10 06:51:19.25 wISdLCFn.net
pa=0.1 ; pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
aa=0:n # 宝くじAを買う候補
Get <- function(a){sum(sample(A0,a))*Wa + sum(sample(B0,n-a))*Wb} #賞金総額
foo <- function(k=10^3){
which.max2 <- function(x){which(x==max(x))-1}
MAX=NULL
MAX.AWD=NULL
for(i in 1:k){
tmp=sapply(aa,Get)
indx=which.max2(tmp)
MAX=c(MAX,indx) # 最大賞金を獲得したときに買ったAの本数の配列
MAX.AWD=c(MAX.AWD,tmp[indx+1]) #最大賞金が複数回あるときは重複して数える
}
indx=which.max2(MAX.AWD)+1
plot(MAX,MAX.AWD,xlim=c(0,n),ylim=c(n*pa*Wa,8*n*pa*Wa),
xlab='A_lots',ylab='Award',col=rgb(.1,.1,.1,.5))
points(MAX[indx],MAX.AWD[indx],pch=19)
luckyA=MAX[indx]
lucky.AWD=MAX.AWD[indx]
cbind(luckyA,lucky.AWD)
}
foo()
145:卵の名無しさん
17/10/10 16:36:50.99 g2TXzqHS.net
> '+'(1,2)
[1] 3
> '*'(2,3)
[1] 6
> a=1:10
> '['(a,7)
[1] 7
'[' って関数なんだな。
146:卵の名無しさん
17/10/10 21:23:23.64 VF4JdnZT.net
>>129
宝くじAは1000本中300本が当たりで当たりは 100万円の賞金、
宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、
とAの期待値＝Bの期待値の3倍のとき
URLﾘﾝｸ(i.imgur.com)
147:卵の名無しさん
17/10/11 06:48:41.73 8dgUXTG5.net
>>139
続きのコードが書き込めないので画像でアップ
URLﾘﾝｸ(imagizer.imageshack.com)
148:卵の名無しさん
17/10/11 07:14:30.46 8dgUXTG5.net
# pa=0.1
pb=0.01 # 当たる確率
Wa=100 ; Wb=1000 # 賞金
n=100 # 買う総本数
foo <- function(pa=0.1,k=10^3){
# k=10^3 # シミュレーション回数
# Aをi買ったときの賞金総額
a <- function(i){sum(rbinom(i, 1, pa) * Wa) + sum(rbinom(100-i, 1, pb) * Wb)}
which.max2 <- function(x){which(x==max(x))-1}
# 最大賞金総額になるAの購入数（複数可）の配列を返す
b <- function(){
b1=NULL
b1=c(b1,which.max2(sapply(0:n, a)))
return(b1)
}
b2=NULL
for(j in 1:k){
b2=c(b2,b())
}
hist(b2,freq=FALSE,breaks=20,col=sample(colors(),1),xlab='',ylab='',yaxt='n',
main=paste('pa =',pa))
lines(density(b2),lwd=1)
print(summary(b2))
print(c(Mode=as.numeric(names(which.max(table(b2)))))) # モード値
print(c(length=length(b2))) # 最大値となる数は複数最大値があるのでｎを超える
invisible(b2)
}
par(mfrow=c(2,2))
for(i in c(0.15,0.2,0.25,0.3)) foo(i,500)
149:卵の名無しさん
17/10/11 20:53:31.06 8dgUXTG5.net
確かにマジックだな
##-- now some "magic" to do the 4 regressions in a loop:
ff <- y ~ x
mods <- setNames(as.list(1:4), paste0("lm", 1:4))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
## or ff[[2]] <- as.name(paste0("y", i))
## ff[[3]] <- as.name(paste0("x", i))
mods[[i]] <- lmi <- lm(ff, data = anscombe)
print(anova(lmi))
}
## See how close they are (numerically!)
sapply(mods, coef)
lapply(mods, function(fm) coef(summary(fm)))
## Now, do what you should have done in the first place: PLOTS
op <- par(mfrow = c(2, 2), mar = 0.1+c(4,4,1,1), oma = c(0, 0, 2, 0))
for(i in 1:4) {
ff[2:3] <- lapply(paste0(c("y","x"), i), as.name)
plot(ff, data = anscombe, col = "red", pch = 21, bg = "orange", cex = 1.2,
xlim = c(3, 19), ylim = c(3, 13))
abline(mods[[i]], col = "blue")
}
mtext("Anscombe's 4 Regression data sets", outer = TRUE, cex = 1.5)
par(op)
150:卵の名無しさん
17/10/11 21:31:48.49 8dgUXTG5.net
[ って関数だったんだ
(d=outer(11:19,11:19,'*'))
d[,3]
apply(d,1,'[',3)
151:卵の名無しさん
17/10/11 23:59:07.55 8dgUXTG5.net
(d=array(1:48, dim=c(6, 4, 2)))
# を，3行2列おきに平均して
# ,,1
# [,1] [,2]
# [1,] 5 17
# [2,] 8 20
# ,,2
# [,1] [,2]
# [1,] 29 41
# [2,] 32 44
#
a=3
b=2
c=1
l=6
m=4
n=2
re=NULL
for(k in 1:(n/c)){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)]))
}
}
}
array(re,dim=c(m/b,l/a,n))
152:卵の名無しさん
17/10/12 18:41:56.42 IVfk5sKl.net
fs<-function(){ #スパゲッティ・ソース
d=array(1:(l*m*n),dim=c(l,m,n))
re=NULL
for(k in 1:(n/c)){
for(j in 1:(m/b)){
for(i in 1:(l/a)){
re=append(re,mean(d[1:a+a*(i-1),1:b+b*(j-1),1:c+c*(k-1)]))
}}}
array(re,dim=c(m/b,l/a,n))
}
fe <- function(){ # Expertソース
A <- array(1:(l*m*n), dim=c(l,m,n))
B <- expand.grid(1:(l %/% a), 1:(m %/% b),1:(n %/% c))
f <- function(i, j, k){mean(A[(a*(i-1)+1):(a*i),(b*(j-1)+1):(b*j),(c*(k-1)+1):(c*k)])}
re=mapply(f, B[,1], B[,2],B[,3])
array(re,dim=c(m%/%b,l%/%a,n%/%c))
}
a=4
b=2
c=1
l=400
m=200
n=2
> system.time(fs())
user system elapsed
2.17 0.00 2.18
> system.time(fe())
user system elapsed
1.14 0.00 1.14
153:卵の名無しさん
17/10/14 04:34:05.51 IFvRV3Gd.net
x=11:19
y=11:19
nx=length(x)
ny=length(y)
re=matrix(rep(NA,nx*ny),ncol=nx)
for(i in 1:nx){
for(j in 1:ny){
re[j,i]=x[i]*y[j]
}
}
re
a=expand.grid(x,y)
b=mapply('*',a[,1],a[,2])
matrix(b,nx)
outer (x,y,'*')
154:卵の名無しさん
17/10/14 06:39:17.90 IFvRV3Gd.net
x=1:999
y=1:999
nx=length(x)
ny=length(y)
fs = function (){
re=matrix(rep(NA,nx*ny),ncol=nx)
for(i in 1:nx){
for(j in 1:ny){
re[j,i]=x[i]*y[j]
}
}
re
}
fe = function (){
a=expand.grid(x,y)
b=mapply('*',a[,1],a[,2])
matrix(b,nx)
}
fo = function (){
outer (x,y,'*')
}
system.time(fs())
system.time(fe())
system.time(fo())
155:卵の名無しさん
17/10/14 19:44:47.68 +G1xzFga.net
##
m=7 ; n=17
Z=matrix(rep(NA,m*n),m,n)
colnames(Z)=paste0('n',0:(n-1))
rownames(Z)=paste0('m',0:(m-1))
f <- function(x) c(x%%m+1,x%%n+1)
for(i in 0:(n*m-1)){
idx=f(i) ; idx
Z[idx[1],idx[2]]=i
}
length(unique(Z))
Z
NZ <- function(x,y,m=7,n=17){
f <- function(x) c(x%%m+1,x%%n+1)
for(i in 0:(n*m-1)){
idx=f(i)
Z[idx[1],idx[2]]=i
}
Z[(x%%m+1),(y%%n+1)]
}
NZ(5,15)
NZ(6,0)
a=expand.grid(0:6,0:16)
res=mapply(NZ,a[,1],a[,2])
Z=matrix(res,m,n)
colnames(Z)=paste0('n',0:(n-1))
rownames(Z)=paste0('m',0:(m-1))
Z
length(unique(Z))
156:卵の名無しさん
17/10/14 19:45:01.23 +G1xzFga.net
where2 <- function(x,m=7,n=17){
modm=paste0('(mod ',m,')')
modn=paste0('(mod ',n,')')
res=c(x%%m,x%%n)
names(res)=c(modm,modn)
return(res)
}
where2(100)
157:卵の名無しさん
17/10/14 19:45:34.01 +G1xzFga.net
## 和で直積
fij <- function(i,j){
ai=where2(i)
aj=where2(j)
aiaj=ai+aj
NZ(aiaj[1],aiaj[2]) == i+j
}
a=expand.grid(0:(m-1),0:(n-1))
mapply(fij,a[,1],a[,2])
## 積で直積
fij <- function(i,j){
ai=where2(i)
aj=where2(j)
aiaj=ai * aj
NZ(aiaj[1],aiaj[2]) == i * j
}
a=expand.grid(0:(m-1),0:(n-1))
mapply(fij,a[,1],a[,2])
158:卵の名無しさん
17/10/14 19:46:23.30 +G1xzFga.net
> mapply(fij,a[,1],a[,2])
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[18] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[35] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[52] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[69] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[86] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
[103] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
同型写像確認
159:卵の名無しさん
17/10/14 22:27:16.20 +G1xzFga.net
NZ <- function(x,y,z, l=13,m=15,n=17){
lmn=l*m*n
N=0:(lmn-1)
.M=rbind(nl=N%%l,nm=N%%m,nn=N%%n)
which(apply((.M - c(x,y,z)),2,function(x) sum(x^2))==0)-1
}
l=13 ; m=15 ; n=17
ll=0:(l-1); mm=0:(m-1) ; nn=0:(n-1)
a=expand.grid(ll,mm,nn)
res=mapply(NZ,a[,1],a[,2],a[,3])
length(unique(res)) == l*m*n
summary(res)
plot(sort(res))
hist(res,breaks=l*m*n)
l*m*nNZ <- function(x,y,z, l=13,m=15,n=17){
lmn=l*m*n
N=0:(lmn-1)
.M=rbind(nl=N%%l,nm=N%%m,nn=N%%n)
which(apply((.M - c(x,y,z)),2,function(x) sum(x^2))==0)-1
}
l=13 ; m=15 ; n=17
ll=0:(l-1); mm=0:(m-1) ; nn=0:(n-1)
a=expand.grid(ll,mm,nn)
res=mapply(NZ,a[,1],a[,2],a[,3])
length(unique(res)) == l*m*n
summary(res)
plot(sort(res))
hist(res,breaks=l*m*n)
l*m*n 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
160:卵の名無しさん
17/10/15 07:15:02.95 hnCh54BK.net
URLﾘﾝｸ(en.wikipedia.org)
を参考にオイラー関数をグラフにしてみた。
URLﾘﾝｸ(i.imgur.com)
phi <- function(n){
nn=0:(n-1)
f=function(x,y) (x*y)%%n
names(nn)=paste0('n',nn)
z=outer(nn,nn,f)
i=which(gmp::gcd(1:(n-1),n)==1)
length(i)
}
#Fourier transform
phi.F <- function(N){
nn=1:N
sum( gmp::gcd(nn,N)*cos(2*pi*nn/N))
}
N=1000
nn=1:N
y=sapply(nn,phi)
plot(nn,y,col=sample(colours()),pch=19,ylab='',xlab='n',main='Euler\'s totient function')
y1=sapply(nn,phi.F)
points(nn,y1,pch='.',cex=1.5)
161:卵の名無しさん
17/10/16 14:46:14.17 SUEY/256.net
(p+1)^(p^(n-1)) ≡1 (mod p^n)
(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))
162:卵の名無しさん
17/10/16 18:24:23.14 SUEY/256.net
# (p+1)^(p^(n-1)) ≡1 (mod p^n)
beki <- function(x,y,p){ # (x^y)%%p
if(y==0) return(1)
if(y==1) return(x%%p)
re=numeric(y)
re[1]=x
for(i in 1:(y-1)) re[i+1] = (x*re[i])%%p
return(re[y])
}
f <- function(p,n){ # (p+1)^(p^(n-1)) ≡1 (mod p^n)
beki(p+1,p^(n-1),p^n)
}
p=(2:9)[gmp::isprime(2:9)!=0]
n=2:9
a=expand.grid(p,n)
> mapply(f,a[,1],a[,2])
> mapply(f,a[,1],a[,2])
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 👀
Rock54: Caution(BBR-MD5:0be15ced7fbdb9fdb4d0ce1929c1b82f)
163:卵の名無しさん
17/10/16 19:48:05.23 SUEY/256.net
(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))
p:奇素数
164:卵の名無しさん
17/10/18 20:51:15.46 r7jfJx+N.net
# (x^y)%%p
beki <- function(x,y,p){
if(!y) return(x^y)
re=numeric()
re[1]=x%%p
for(i in 1:y) re[i+1] = (x*re[i])%%p
return(re[y])
}
# ｎのmod p での位数を求める
isu <- function(n,p){
re=numeric()
for(i in 1:(p-1)){
re[i]=beki(n,i,p)
}
which.min(re)
}
165:卵の名無しさん
17/10/19 15:04:12.35 lylsDJ3i.net
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
e =c(1,2,3)
d =c(3,1,2)
d2=c(2,3,1)
t =c(1,3,2)
td=c(3,2,1)
td2=c(2,1,3)
(.M=rbind(e,d,d2,t,td,td2))
(names=rownames(.M))
equal <- function(x,y) sum((x-y)^2)==0
f <- function(x){ #
re=NULL
for(i in 1:nrow(.M)){
if(equal(.M[i,],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
166:卵の名無しさん
17/10/19 15:04:48.68 lylsDJ3i.net
.L = list(e,d,d2,t,td,td2)
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),nrow(.M)))
.M1=matrix(names[.a], nrow(.M))
rownames(.M1)=names
colnames(.M1)=names
> print(.M1, quote=FALSE)
e d d2 t td td2
e e d d2 t td td2
d d d2 e td2 t td
d2 d2 e d td td2 t
t t td td2 e d d2
td td td2 t d2 e d
td2 td2 t td d d2 e
167:卵の名無しさん
17/10/19 16:00:22.26 lylsDJ3i.net
結合側の検証
k <- function(x,y,z) equal(tikan2(x,tikan2(y,z)), tikan2(tikan2(x,y),z))
(a=expand.grid(.L,.L,.L))
b=mapply(k,a[,1],a[,2],a[,3])
> summary(b)
Mode TRUE
logical 216
168:卵の名無しさん
17/10/19 19:20:37.10 pl4P+6GE.net
>>162
結合則
169:卵の名無しさん
17/10/19 19:26:48.29 ZZphLwM2.net
演算が閉じている、単位元が存在する、逆元が存在する、の確認は容易だけど
結合則の確認は手計算だと手間がかかる。
プログラムを組んで確認できて納得。
170:卵の名無しさん
17/10/20 13:57:10.48 hxLmqOB4.net
紛らわしいなぁ
群論における「位数」とは、二つの異なる定義があり、それぞれ、
①有限群Gの位数＝有限群Gの元の個数、
②有限群Gの元aの位数＝有限群Gの元aにおいて、a^m=e となる最小の正の整数、
となります。
171:卵の名無しさん
17/10/20 14:07:59.55 hxLmqOB4.net
## 位数２ｎの二面体群の演算表を作る
# D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
n=7
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
e=1:n
d1=c(n,1:(n-1)) # 360°/n 回転移動
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i in 1:(n-1)){
Mnd[i+1,]=tikan(Mnd[i,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
172:卵の名無しさん
17/10/20 14:08:06.73 hxLmqOB4.net
t=n:1 ## 対称移動
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i in 1:(n-1)){
Mnt[i+1,]=tikan(Mnt[i,],d1)
}
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mnt
(Mn=rbind(Mnd,Mnt))
names=rownames(Mn)
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i in 1:(2*n)) .L[[i]]=Mn[i,]
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),length(.L)))
.M1=matrix(names[.a], length(.L))
rownames(.M1)=names
colnames(.M1)=names
.M1
print(.M1, quote=FALSE)
173:卵の名無しさん
17/10/20 14:08:34.26 hxLmqOB4.net
> print(.M1, quote=FALSE)
e d1 d2 d3 d4 d5 d6 t td1 td2 td3 td4 td5 td6
e e d1 d2 d3 d4 d5 d6 t td1 td2 td3 td4 td5 td6
d1 d1 d2 d3 d4 d5 d6 e td6 t td1 td2 td3 td4 td5
d2 d2 d3 d4 d5 d6 e d1 td5 td6 t td1 td2 td3 td4
d3 d3 d4 d5 d6 e d1 d2 td4 td5 td6 t td1 td2 td3
d4 d4 d5 d6 e d1 d2 d3 td3 td4 td5 td6 t td1 td2
d5 d5 d6 e d1 d2 d3 d4 td2 td3 td4 td5 td6 t td1
d6 d6 e d1 d2 d3 d4 d5 td1 td2 td3 td4 td5 td6 t
t t td1 td2 td3 td4 td5 td6 e d1 d2 d3 d4 d5 d6
td1 td1 td2 td3 td4 td5 td6 t d6 e d1 d2 d3 d4 d5
td2 td2 td3 td4 td5 td6 t td1 d5 d6 e d1 d2 d3 d4
td3 td3 td4 td5 td6 t td1 td2 d4 d5 d6 e d1 d2 d3
td4 td4 td5 td6 t td1 td2 td3 d3 d4 d5 d6 e d1 d2
td5 td5 td6 t td1 td2 td3 td4 d2 d3 d4 d5 d6 e d1
td6 td6 t td1 td2 td3 td4 td5 d1 d2 d3 d4 d5 d6 e
174:卵の名無しさん
17/10/20 14:21:16.31 hxLmqOB4.net
これが抜けていた。
f <- function(x){
re=NULL
for(i in 1:length(.L)){
if(equal(.L[[i]],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
175:卵の名無しさん
17/10/20 14:30:53.07 hxLmqOB4.net
## 位数２ｎの二面体群の演算表を作る
# D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
n=7
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
equal <- function(x,y) sum((x-y)^2)==0
176:卵の名無しさん
17/10/20 14:32:25.05 hxLmqOB4.net
f <- function(x){
re=NULL
for(i in 1:length(.L)){
if(equal(.L[[i]],x)) return(i) # i for names[i]
}
}
g <- function(x,y){
f(tikan2(x,y)) # first y, then x makes names[i]
}
e=1:n
d1=c(n,1:(n-1)) # 360°/n 回転移動
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i in 1:(n-1)){
Mnd[i+1,]=tikan(Mnd[i,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
177:卵の名無しさん
17/10/20 14:32:49.69 hxLmqOB4.net
t=n:1 ## 対称移動
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i in 1:(n-1)){
Mnt[i+1,]=tikan(Mnt[i,],d1)
}
rownames(Mnt)=c('t',paste0('td',1:(n-1)))
Mnt
(Mn=rbind(Mnd,Mnt))
names=rownames(Mn)
.L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1)
for(i in 1:(2*n)) .L[[i]]=Mn[i,]
names(.L)=names
.L
a=expand.grid(.L,.L)
(.a=mapply(g,a[,1],a[,2]))
(.M0=matrix(mapply(g,a[,1],a[,2]),length(.L)))
.M1=matrix(names[.a], length(.L))
rownames(.M1)=names
colnames(.M1)=names
.M1
print(.M1, quote=FALSE)
178:卵の名無しさん
17/10/20 14:41:11.07 hxLmqOB4.net
>170-172で動作確認
179:卵の名無しさん
17/10/21 07:30:34.81 H3fe+vIj.net
n=10での演算表
e d1 d2 d3 d4 d5 d6 d7 d8 d9 t td1 td2 td3 td4 td5 td6 td7 td8 td9
e e d1 d2 d3 d4 d5 d6 d7 d8 d9 t td1 td2 td3 td4 td5 td6 td7 td8 td9
d1 d1 d2 d3 d4 d5 d6 d7 d8 d9 e td9 t td1 td2 td3 td4 td5 td6 td7 td8
d2 d2 d3 d4 d5 d6 d7 d8 d9 e d1 td8 td9 t td1 td2 td3 td4 td5 td6 td7
d3 d3 d4 d5 d6 d7 d8 d9 e d1 d2 td7 td8 td9 t td1 td2 td3 td4 td5 td6
d4 d4 d5 d6 d7 d8 d9 e d1 d2 d3 td6 td7 td8 td9 t td1 td2 td3 td4 td5
d5 d5 d6 d7 d8 d9 e d1 d2 d3 d4 td5 td6 td7 td8 td9 t td1 td2 td3 td4
d6 d6 d7 d8 d9 e d1 d2 d3 d4 d5 td4 td5 td6 td7 td8 td9 t td1 td2 td3
d7 d7 d8 d9 e d1 d2 d3 d4 d5 d6 td3 td4 td5 td6 td7 td8 td9 t td1 td2
d8 d8 d9 e d1 d2 d3 d4 d5 d6 d7 td2 td3 td4 td5 td6 td7 td8 td9 t td1
d9 d9 e d1 d2 d3 d4 d5 d6 d7 d8 td1 td2 td3 td4 td5 td6 td7 td8 td9 t
t t td1 td2 td3 td4 td5 td6 td7 td8 td9 e d1 d2 d3 d4 d5 d6 d7 d8 d9
td1 td1 td2 td3 td4 td5 td6 td7 td8 td9 t d9 e d1 d2 d3 d4 d5 d6 d7 d8
td2 td2 td3 td4 td5 td6 td7 td8 td9 t td1 d8 d9 e d1 d2 d3 d4 d5 d6 d7
td3 td3 td4 td5 td6 td7 td8 td9 t td1 td2 d7 d8 d9 e d1 d2 d3 d4 d5 d6
td4 td4 td5 td6 td7 td8 td9 t td1 td2 td3 d6 d7 d8 d9 e d1 d2 d3 d4 d5
td5 td5 td6 td7 td8 td9 t td1 td2 td3 td4 d5 d6 d7 d8 d9 e d1 d2 d3 d4
td6 td6 td7 td8 td9 t td1 td2 td3 td4 td5 d4 d5 d6 d7 d8 d9 e d1 d2 d3
td7 td7 td8 td9 t td1 td2 td3 td4 td5 td6 d3 d4 d5 d6 d7 d8 d9 e d1 d2
td8 td8 td9 t td1 td2 td3 td4 td5 td6 td7 d2 d3 d4 d5 d6 d7 d8 d9 e d1
td9 td9 t td1 td2 td3 td4 td5 td6 td7 td8 d1 d2 d3 d4 d5 d6 d7 d8 d9 e
180:卵の名無しさん
17/10/21 19:00:23.36 TcIT7+c6.net
## D6で位数3の部分群を虱潰しに探す
a=t(combn(1:12,3))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
# names(H)=names(G[c(a[ii,1],a[ii,2],a[ii,3])])
M=matrix(NA,length(H),length(G))
for(i in 1:length(H)){
for(j in 1:length(G)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
#colnames(M)=paste0(names(G),'*')
#rownames(M)=names(H)
#print(M,quote=FALSE)
re=matrix(NA,nrow(M),ncol(M))
for(i1 in 1:ncol(M)){
re[,i1]=sort(M[,i1])
}
#print(re,quote=FALSE)
#print(unique(t(re)),quote=FALSE)
RES[ii] =nrow(unique(t(re)))
}
RES
idx=which(RES==4)
(b=a[idx,])
print(matrix(names[b],ncol=3),quote=FALSE)
181:卵の名無しさん
17/10/21 19:00:49.55 TcIT7+c6.net
> (b=a[idx,])
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
[3,] 7 9 11
[4,] 8 10 12
> print(matrix(names[b],ncol=3),quote=FALSE)
[,1] [,2] [,3]
[1,] e d2 d4
[2,] d1 d3 d5
[3,] t td2 td4
[4,] td1 td3 td5
182:卵の名無しさん
17/10/21 20:08:08.06 TcIT7+c6.net
## 二面体群Dnで位数mの部分群を虱潰しに探す
n=9 ; m=3
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f1 <- function(x,L){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(names(L)[i]) # i for names[i]
}
}
183:卵の名無しさん
17/10/21 20:19:32.79 TcIT7+c6.net
URLﾘﾝｸ(i.imgur.com)
184:卵の名無しさん
17/10/21 20:20:03.91 TcIT7+c6.net
##
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
idx=which(RES==(2*n/m)) # Lagrangea's theorem
(b=a[idx,])
print(matrix(names[b],ncol=m),quote=FALSE)
185:卵の名無しさん
17/10/21 20:21:10.03 TcIT7+c6.net
> (b=a[idx,])
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
[4,] 10 13 16
[5,] 11 14 17
[6,] 12 15 18
> print(matrix(names[b],ncol=m),quote=FALSE)
[,1] [,2] [,3]
[1,] e d3 d6
[2,] d1 d4 d7
[3,] d2 d5 d8
[4,] t td3 td6
[5,] td1 td4 td7
[6,] td2 td5 td8
# 御明算
186:卵の名無しさん
17/10/21 21:11:33.21 TcIT7+c6.net
>>179
(debugged)
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
idx=which(RES==(2*n/m)) # Lagrangea's theorem
b=a[idx,]
(B=b[b[,1]==1,])
print(matrix(names[B],ncol=m),quote=FALSE)
187:卵の名無しさん
17/10/21 23:55:34.83 TcIT7+c6.net
## 二面体群Dnですべての部分群を返す
dihedral_group <- function(n=9,m=3){ # start
tikan <- function(X,Y=.Y){ # e.g. c(1,2,3) -> .Y=c(2,1,3)
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f1 <- function(x,L){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(names(L)[i]) # i for names[i]
}
}
g12 <- function(x,y,L=G){
f1(tikan2(x,y),L)
}
equal <- function(x,y) sum((x-y)^2)==0
188:卵の名無しさん
17/10/21 23:56:01.41 TcIT7+c6.net
e=1:n
d1=c(n,1:(n-1)) # rotation
Mnd=matrix(NA,n,n)
Mnd[1,]=e
for(i0 in 1:(n-1)){
Mnd[i0+1,]=tikan(Mnd[i0,],d1)
}
rownames(Mnd)=c('e',paste0('d',1:(n-1)))
Mnd
t=n:1 ## mirror
tikan(t,d1) # D2n( e,d1,d2,d3,..,dn-1, t,td,td2,td3,...,tdn-1)
Mnt=matrix(NA,n,n)
Mnt[1,]=t
for(i2 in 1:(n-1)){
Mnt[i2+1,]=tikan(Mnt[i2,],d1)
}
189:Mnt rownames(Mnt)=c('t',paste0('td',1:(n-1))) Mn=rbind(Mnd,Mnt) names=rownames(Mn) Mn .L=list(2*n) # list(Mn[1,],Mn[2,],...Mn[,n]) list(e,d1,d2,..,dn-1,t,td1,td2,..,tdn-1) for(i3 in 1:(2*n)) .L[[i3]]=Mn[i3,] names(.L)=names G=.L
190:卵の名無しさん
17/10/21 23:56:34.66 TcIT7+c6.net
a=t(combn(1:(2*n),m))
RES=numeric(nrow(a))
for(ii in 1:nrow(a)){
H=list() # H=list(G[[a[ii,1]]],G[[a[ii,2]]],G[[a[ii,3]]])
for(jj in 1:m){
H[[jj]]=G[[a[ii,jj]]]
}
M=matrix(NA,m,2*n)
for(i in 1:m){
for(j in 1:(2*n)){
M[i,j] = g12(G[[j]],H[[i]])
}
}
re=matrix(NA,m,2*n)
for(i1 in 1:(2*n)){
re[,i1]=sort(M[,i1])
}
RES[ii] =nrow(unique(t(re)))
}
RES
idx=which(RES==(2*n/m)) # Lagrangea's theorem
(b=a[idx,])
#b=as.matrix(b)
if(m==1||m==2*n){
B=1:m
}else{
B=b[b[,1]==1,]
}
matrix(names[B],ncol=m)
} # End of Function, dihedral_group
191:卵の名無しさん
17/10/21 23:57:03.96 TcIT7+c6.net
.print <- function(x) print(x,quote=FALSE)
sub_group <- function(n){
N=2*n
xx=which(N%%(1:N)==0)
res=lapply(xx,function(x)dihedral_group(n,x))
names(res)=xx
.print(res)
}
192:卵の名無しさん
17/10/21 23:57:45.77 TcIT7+c6.net
> sub_group(6)
$`1`
[,1]
[1,] e
$`2`
[,1] [,2]
[1,] e d3
[2,] e t
[3,] e td1
[4,] e td2
[5,] e td3
[6,] e td4
[7,] e td5
$`3`
[,1] [,2] [,3]
[1,] e d2 d4
$`4`
[,1] [,2] [,3] [,4]
[1,] e d3 t td3
[2,] e d3 td1 td4
[3,] e d3 td2 td5
$`6`
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] e d1 d2 d3 d4 d5
[2,] e d2 d4 t td2 td4
[3,] e d2 d4 td1 td3 td5
$`12`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] e d1 d2 d3 d4 d5 t td1 td2 td3 td4 td5
193:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 09:18:14.91 bCTtjSj1.net
# Greatest Common Divisor with Euclidean algorithm
GCD <- function(x,y){
if(round(x)!=x | round(y)!=y) stop('Not integer')
if(!x&!y) return(NA)
a=max(c(abs(x),abs(y)))
b=min(c(abs(x),abs(y)))
if(!a*b) return(a)
r=integer()
r[1]=a
r[2]=b
i=1
r[i+2]=r[i]%%r[i+1]
while(r[i+2]!=0 & r[i+2]!=1){
i=i+1
r[i+2]=r[i]%%r[i+1]
}
return(ifelse(r[i+2]==0,r[i+1],1))
}
GCD(2.3,4)
GCD(0,0)
GCD(24,15)
GCD(16,15)
GCD(5,0)
GCD(24,-15)
GCD(-24,-15)
194:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 13:48:21.69 bCTtjSj1.net
> sub_group(9)
$`1`
[,1]
[1,] e
$`2`
[,1] [,2]
[1,] e t
[2,] e td1
[3,] e td2
[4,] e td3
[5,] e td4
[6,] e td5
[7,] e td6
[8,] e td7
[9,] e td8
$`3`
[,1] [,2] [,3]
[1,] e d3 d6
$`6`
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] e d3 d6 t td3 td6
[2,] e d3 d6 td1 td4 td7
[3,] e d3 d6 td2 td5 td8
$`9`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] e d1 d2 d3 d4 d5 d6 d7 d8
$`18`
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15]
[1,] e d1 d2 d3 d4 d5 d6 d7 d8 t td1 td2 td3 td4 td5
[,16] [,17] [,18]
[1,] td6 td7 td8
195:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 15:51:11.54 bCTtjSj1.net
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
# Gを有限群とすると g∈G n=|G| （位数） g^n=e :単位元
n=length(G)
re=list()
for(i in 1:n){
tikan2(G[[i]],G[[i]])
for(j in 1:n) data[[j]]=G[[i]]
re[[i]]=Reduce(tikan2,data) # Reduce(function(x,y)x+y,c(1,2,3,4),accumulate = TRUE)
}
> re
[[1]]
[1] 1 2 3 4 5 6 7 8 9
[[2]]
[1] 1 2 3 4 5 6 7 8 9
......
[[17]]
[1] 1 2 3 4 5 6 7 8 9
[[18]]
[1] 1 2 3 4 5 6 7 8 9
196:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 16:46:17.03 bCTtjSj1.net
> Reduce(function(x,y)x+y,c(1,2,3,4,5),accumulate = TRUE)
[1] 1 3 6 10 15
> Reduce('+',c(1,2,3,4,5),accu=TRUE) ; cumsum(1:5)
[1] 1 3 6 10 15
[1] 1 3 6 10 15
> Reduce('*',c(1,2,3,4,5)) ; factorial(5)
[1] 120
[1] 120
> beki2 <- function(x,y,p,...){ #x^y%%p
+ if(y==0) return(1%%p)
+ data=rep(x,y)
+ Reduce(function(x,y) (x*y)%%p, data,...)
+ }
> beki2(2,10,10,accumulate=TRUE)
[1] 2 4 8 6 2 4 8 6 2 4
> beki2(2,100,10)
[1] 6
> beki2(2,0,10)
[1] 1
197:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 17:29:41.55 bCTtjSj1.net
# When p is a prime number and a is coprime to p,
#
# a^(p-1)≡1 (mod p)
#
# Check if it works when p is lower than N.
N=100
(p100=(1:N)[!!gmp::isprime(1:N)] ) # primes < 100
m=length(p100)
cop=list()
for(i in 1:m){
p=p100[i]
cop[[i]]=(1:p)[gmp::gcd(1:p,p)==1] # disjoint,coprime to p100[i]
}
names(cop)=p100
cop
re=cop
for(i in 1:m){
p=p100[i]
for(j in 1:length(cop[[i]])){
re[[i]][j] = beki2(cop[[i]][j],p-1,p) # coprime^(p-1)%%p
}
}
re
> str(re)
List of 25
$ 2 : int 1
$ 3 : int [1:2] 1 1
...
$ 89: int [1:88] 1 1 1 1 1 1 1 1 1 1 ...
$ 97: int [1:96] 1 1 1 1 1 1 1 1 1 1 ...
フェルマーの定理
198:名無しさん＠そうだ選挙に行こう！ Go to vote!
17/10/22 18:21:14.08 r/LJPTCR.net
>>191
フェルマーも�
199:ｬ定理の方ね
200:卵の名無しさん
17/10/25 19:11:29.66 nTDNW+TM.net
>>122
##
n=10
conf.level=0.95
PMF <- function(x) x^n # AUC==integrate(PMF,0,1) == 1/(n+1)
PDF <- function(x) (n+1)*x^n # PDF==PMF/AUC
Ex= (n+1)/(n+2) # integrate(x*PDF,0,1)==integrate((n+1)*x^(n+1),0,1)== (n+1)/(n+2)
CDF <- function(x) x^(n+1) # == integrate(PDF,0,x)
lwr = (1-conf.level)^(1/(n+1)) # CDF(lwr) == 1-conf.level
exp(log(1-conf.level)/(n+1))
curve(PDF)
curve(x*PDF(x))
integrate(function(x) x*PDF(x),0,1)$value
integrate(function(x) x*(n+1)*x^n,0,1)$value # integrate(function(x)(n+1)*x^(n+1),0,1)
(n+1)/(n+2)
# |((n+1)/(n+2))*x^(n+2)|[0,1]
201:卵の名無しさん
17/10/25 19:13:22.53 nTDNW+TM.net
>>193
## binomとの比較
binom::binom.confint(c(1,10,100),c(1,10,100),conf=0.95)
# n out of n
# mean=(n+1)/(n+2)
# lower=(n+1)√(1-conf.level) = 0.05^(1/(n+1))
non=function(n,conf.level=0.95){ # n hits out of n shoots
b=binom::binom.confint(n,n)[,c('method','mean','lower')]
n_1_root=function(n)(0.05)^(1/(n+1))
a=data.frame(method='root(n+1)',mean=(n+1)/(n+2),lower=(1-conf.level)^(1/(n+1)))
rbind(a,b)
}
> non(1)
method mean lower
1 root(n+1) 0.6666667 0.22360680
2 agresti-coull 1.0000000 0.16749949
3 asymptotic 1.0000000 1.00000000
4 bayes 0.7500000 0.22851981
5 cloglog 1.0000000 0.02500000
6 exact 1.0000000 0.02500000
7 logit 1.0000000 0.02500000
8 probit 1.0000000 0.02500000
9 profile 1.0000000 0.13811125
10 lrt 1.0000000 0.14650325
11 prop.test 1.0000000 0.05462076
12 wilson 1.0000000 0.20654931
202:卵の名無しさん
17/10/25 19:13:48.14 nTDNW+TM.net
> non(10)
method mean lower
1 root(n+1) 0.9166667 0.7615958
2 agresti-coull 1.0000000 0.6791127
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9545455 0.8292269
5 cloglog 1.0000000 0.6915029
6 exact 1.0000000 0.6915029
7 logit 1.0000000 0.6915029
8 probit 1.0000000 0.6915029
9 profile 1.0000000 0.7303058
10 lrt 1.0000000 0.8252466
11 prop.test 1.0000000 0.6554628
12 wilson 1.0000000 0.7224672
> non(100)
method mean lower
1 root(n+1) 0.9901961 0.9707748
2 agresti-coull 1.0000000 0.9555879
3 asymptotic 1.0000000 1.0000000
4 bayes 0.9950495 0.9810231
5 cloglog 1.0000000 0.9637833
6 exact 1.0000000 0.9637833
7 logit 1.0000000 0.9637833
8 probit 1.0000000 0.9637833
9 profile 1.0000000 0.9670434
10 lrt 1.0000000 0.9809757
11 prop.test 1.0000000 0.9538987
12 wilson 1.0000000 0.9630065
203:卵の名無しさん
17/10/25 19:24:32.39 nTDNW+TM.net
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
URLﾘﾝｸ(i.imgur.com)
204:卵の名無しさん
17/10/25 19:24:57.43 nTDNW+TM.net
>>195
n=10 の時のPDFと９５％CI下限値0.7615958
URLﾘﾝｸ(i.imgur.com)
205:卵の名無しさん
17/10/25 20:40:32.01 nTDNW+TM.net
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}
206:卵の名無しさん
17/10/26 12:27:43.35 ljbBi1cA.net
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.auｎc(1)-f.auc(0)) ; 1/auc
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
207:卵の名無しさん
17/10/26 16:21:26.13 ljbBi1cA.net
# ｆ回失敗した後にｓ回目の成功,K:シミュレーション回数
Dotsubo <- function(s=1,f=4,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(p) choose(s+f-1,f)*(1-p)^f*p^s
AUC= (gamma(s+f)*gamma(s+1)) / (gamma(s)*gamma(f+s+2))
PDF <- function(p) (1-p)^f*p^s*gamma(f+s+2)/(gamma(f+1)*gamma(s+1))
curve(PDF)
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),2),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
208:卵の名無しさん
17/10/26 17:43:49.64 z4lFiDnk.net
ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
209:卵の名無しさん
17/10/27 15:54:12.90 dzxKDqmi.net
# ｎ発ｒ中の狙撃手がＮ発狙撃するときの命中数を返す
Golgo.sim <- function(N, n, r, k=10^3,Print=TRUE){ # k:シミュレーション回数
f <-function(S,.N=N,.n=n){ # 成績サンプル：命中S個、外れ(N-S)個
y=c(rep(1,S),rep(0,.N-S))
sum(sample(y,.n)) # その成績サンプルからｎ個数取り出したときの命中数
}
xx=r:N # ｒ未満ではｒ個命中することはないので除外
SS=NULL # 容れ子
for(i in 1:k){
x=sapply(xx,f) # 命中数の配列
SS=append(SS,which(x==r)-1+r) # 命中数がｒ個のときの成績サンプルの命中数Ｓの配列をつくる
}
print(summary(SS))
print(quantile(SS,probs=c(0.025,0.05,0.95,0.975)))
print(c(mode=names(which.max(table(SS)))),quote=FALSE)
if(Print) {
hist(SS,xlim=c(0,N),freq=FALSE,col=sample(colors(),1),main='',xlab='Hits')
lines(density(SS))}
invisible(SS)
}
210:卵の名無しさん
17/10/27 15:55:45.80 pwjeiI6l.net
# n発r中の期待値
Golgo <- function(n=3,r=1,cl=0.95,K=10^6,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
AUC=integrate(PMF,0,1)$value
# library(hypergeo)
# f.auc <- function(x) choose(n,r)*x^(r+1)/(r+1)*hypergeo(r+1,r-n,r+2,x)
# auc=as.numeric(f.auc(1)-f.auc(0))
PDF <- function(x)PMF(x)/AUC
Ex=integrate(function(x)x*PDF(x),0,1)$value
mode=optimize(PDF,c(0,1),maximum = TRUE)$maximum
CDF <- function(x) integrate(PDF,0,x)$value
CDFu0 <- function(x,u0=.025) CDF(x)-u0
lwr=uniroot(CDFu0,u0=(1-cl)/2,c(0,1))$root
upr=uniroot(CDFu0,u0=1-(1-cl)/2,c(0,1))$root
print(c(lower=lwr,mode=mode,mean=Ex,upper=upr))
xx=seq(0,1,by=0.0001)
pp=sapply(xx,PDF)
yy=sample(xx,K,replace=TRUE,prob=pp)
if(Print){hist(yy,freq=FALSE,main='',xlab='Probability',
col=sample(colors(),1),breaks=30)
lines(density(yy),col='gray',lty=3)
curve(PDF,add=TRUE)}
print(quantile(yy,probs=c(.025,.05,.5,.95,.975)))
myy=mean(yy)
dens <- density(yy)
mdyy=dens$x[which.max(dens$y)]
print(c(mode=mdyy,mean=myy))
invisible(yy)
}
211:卵の名無しさん
17/10/29 17:45:50.97 LlbU36d2.net
data{// coin5.stan
int N;
int<lower=0,upper=1> Y[N];
}
parameters{
real<lower=0,upper=1> p;
}
model{
for(n in 1:N)
Y[n] ~ bernoulli(p);
}
data{// coin1.stan
int<lower=0,upper=1> Y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
Y ~ bernoulli(p);
}
212:卵の名無しさん
17/10/29 22:16:10.95 LlbU36d2.net
# クラインの4元群Vが正6面体群S(P6)の正規部分群であることの確認
equal <- function(x,y) sum((x-y)^2)==0
tikan2 <- function(X,Y){ # replace Y first then X
n=length(X)
Z=rep(NA,n)
I=1:n
for(i in 1:n){
j <- which(I==X[i])
Z[i] <- Y[j]
}
return(Z)
}
f2 <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
a =c(2,1,4,3)
b =c(3,4,1,2)
c =c(4,3,2,1)
e =c(1,2,3,4)
d =c(1,3,4,2)
t =c(1,2,4,3)
d2 =tikan2(d,d) # 1 4 2 3 # d3=tikan2(d2,d) == e
td =tikan2(t,d) # 1 3 2 4
td2=tikan2(td,d)
213:卵の名無しさん
17/10/29 22:17:23.29 LlbU36d2.net
V=list(e,a,b,c)
D3=list(e,d,d2,t,td,td2)
gr=expand.grid(V,D3)
.m=mapply(tikan2,gr[,1],gr[,2])
t.m=t(.m)
G=list()
for(i in 1:nrow(t.m)) G[[i]]=t.m[i,]
names(G)=paste0('t',1:length(G))
G # G: S(P6)
lG=length(G)
V
lV=length(V)
.M1=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M1[j,i]=f2(tikan2(G[[i]],V[[j]]))
}
}
.M2=matrix(NA,nrow=lV,ncol=lG)
for (i in 1:lG){
for(j in 1:lV){
.M2[j,i]=f2(tikan2(V[[j]],G[[i]]))
}
}
identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
> identical(apply(.M1,2,sort) ,apply(.M2,2,sort))
[1] TRUE