臨床統計もおもしろいですよ、その１ [無断転載禁止]©2ch.net

**卵の名無しさん** · 2017/05/03(水) 20:04:54.62

　
　内科認定医受験の最低限の知識、
　製薬会社の示してくる臨床データ、
　論文の考察、
　論文を書くときの正当性、
　というのが、臨床統計の今までの目的の大きい部分でしたが、
　
　AI=機械学習の基本も、結局は統計学と確率に支配されます。
　そういう雑多な話をするスレです。
　

**卵の名無しさん** · 2017/05/03(水) 20:09:27.30

　
　明日、正午から1時間お時間のある先生には、BSで放送大学の心理統計学の第二回を
　見ていただくと、モンティ・ホール問題という常識をひっくり返す話があると思います。
　これで、私も混乱しているので、一からやり直そうかと。
　
　区間推計って1/√n-1で決まっていたんじゃないのかな、いえ、1/√nでもいいんじゃないか、
　それだけって私にはお手上げなので...お願いします。
　

**卵の名無しさん** · 2017/05/03(水) 20:13:44.08

　
　もちろん、いつも見ている2x2のΧ^2分析だけ話していてもおもしろくはないでしょ。
　臨床医の知的トレーニング、サイコロでもRの話、何でも雑談してください。
　

**卵の名無しさん** · 2017/05/03(水) 23:40:48.38

最悪の精神科医　古根高 [転載禁止](c)2ch.net
http://potato.2ch.net/test/read.cgi/hosp/1439931587/

過去ログだがブラウザで読める

**卵の名無しさん** · 2017/05/04(木) 10:37:21.49

スレ立て乙。

区間推定っていろんな方法があるね。
正規分布近似がよく用いられるけど。
MotulskyのIntuitive Biostastisticsには正規分布でなくｔ分布近似しろ、
正規分布近似のEXCELは間違いだという趣旨のことが書いてある。
ｔ分布の自由度を∞にすれば正規分布なのでさほど気にすることもないと思っている。
同様に√nも√n-1も大差がないと思っている。
それが結論を左右するならnが小さすぎだと思う。

Rのパッケージbinomに
binom.confint という関数があって種々の信頼区間を吐き出してくれて便利。極端な値でなければどれも一致する。

**卵の名無しさん** · 2017/05/04(木) 15:10:54.15

>>5
　私には、どういう分布を採用していいのかという理論のところが分からなくて、
　Excelなり、Rなりに後はデータを入れてしまう、認定試験も2x2の表が
　まあまあ理解できないと分からないレベルでして（笑）。
　
　正規分布を使うべきか、t分布を使うべきか、どうやって決められるのですか？
　後、不偏分散と分散の差も分かってません（苦笑）。
　

**卵の名無しさん** · 2017/05/04(木) 15:13:12.06

>>5
　ああ、二項分布やポアソン分布だったか、n → ∞になれば正規分布になるという理屈は
　分かるのですが。
　

**卵の名無しさん** · 2017/05/04(木) 15:18:16.74

>>6の続き、
　空腹で低血糖のため、日本語がおかしいですが、
　　2x2の有意差検定の表がまあまあ理解できるレベルで...
　という意味です。
　
　大学の教養学部の時代に、もう30年前の話ですが、分布関数を考えるときには、
　積率母関数を考えなさいと教わりました。たぶん、ラプラス変換と同じで実数関数の
　世界で収めたいという暖かい心遣いの講師の先生だったと思うのですが、今なら
　虚数関数の特性関数の方が使いやすいなと思っております。基礎は堅いのですが
　応用をどうして良いのか、迷ってます。
　

**卵の名無しさん** · 2017/05/04(木) 16:02:02.70

　
　土曜日に5時まで医師としての仕事をしたとします。
　あと1時間の先生も多いですよね。

　で、他方、東京の普通の本社勤務のクリエイティブな仕事をしている人は、9連休でも経済的には
　OKなんですよね？
　
　医師はQOLを論じて議論するわけでしょ、統計学的に明らかに騙されていると思いますが。
　

**卵の名無しさん** · 2017/05/04(木) 17:24:09.15

>>7
ｎ=∞のtの分位関数は正規分布の分位関数と同じという話を書いたつもり。

Rで書けば
> qnorm(0.975)
[1] 1.959964
> qt(0.975,Inf)
[1] 1.959964

というだけの話。
ｔ分布を使う方が区間の幅が広くなるけど。

**卵の名無しさん** · 2017/05/04(木) 18:10:15.18

>>6
どの分布にしたがうかは経験則ではないかと思ってます。

たとえば所得の分布が対数正規分布に従うというのは理論的帰結ではないのだろうと。

過分散を避けるためにポアソン分布でなく負の二項分布を使うとか読んだけど理解できていない。

それが、わからぬままに

基礎からのベイズ統計学の　続編

実践ベイズモデリング -解析技法と認知モデル

も読んだけど

さまざまなモデルが載っていたけれどモデルからの推論が現実と合致していればそれでいいのだろうと深く追求しておりません。

κ値（検査間や検者間の一致の指標）のMCMCサンプルの方法

和が1になる変数simplexと多項分布を使うことを思いつかなかったから、stanのスクリプトをみて眼から鱗。

これでκ値の信頼区間がだせる。

ちなみにMcGee　Evidence-Based Physical Diagnosisの新版を買うかどうか迷っている。

旧版は持っているのでどれほど中身が変わったかな。　Google Bookで立ち読みしたけど髄膜炎のjolt accentuationのκ値の記載はなかったな。

**卵の名無しさん** · 2017/05/04(木) 18:52:09.21

>>8
>2x2の有意差検定のMCMCでのやり方が

基礎からのベイズ統計学の8章に書いてあったので

RのパッケージEpiのtwoby2とか、fisher.testとかprop.testが吐き出す信頼区間と一致するかで遊んだことがある。

シンプソンのパラドックスによく引用される男女の合格率の差

https://www.refsmmat.com/posts/2016-05-08-simpsons-paradox-berkeley.html

を使ってやってみたのが

http://egg.2ch.net/test/read.cgi/hosp/1428282054/759-761

今の俺の知識ではRのパッケージでできることを無理やりMCMCでやっている気がする。

まあ、サンプルからパラメータの信頼区間を求められるのが長所かな。

**卵の名無しさん** · 2017/05/05(金) 01:42:46.64

>>11
　逆にきれいに従うモデルがあるとすれば、そこにこそ真実があるのかも知れませんね。
　遅くなっちゃったので今日はこれまでと言うことで。
　

**卵の名無しさん** · 2017/05/05(金) 11:45:29.28

>>12
　Rに慣れたらやってみます。パッケージすら、選び方がなかなかで...。
　ベイズ統計学の講義を視聴しているのですが、問題は事前分布が分からないような場合に
　どうするのかと言うことが、MCMC法を使う理由のようですね。
　線虫とかイヌの嗅覚で、がんを見つけるときに、母集団の有病率の事前分布をさてどうやって
　規定すればいいのやら、ってのことではないのでしょうか。
　
　朝から木を切っていて疲れました、これは剛体の物理学通りに落ちてくれるのですが、
　予期しない落ち方もありますね（笑）。
　

**卵の名無しさん** · 2017/05/05(金) 11:49:19.71

>>11
　英語で考えると日本語の3倍ぐらい時間がかかってしまうようになってしまって、
　カルテも電子カルテで情報共有と言われると、日本語ばかりになってしまうわけで
　英語で考えたり、読んだりが大変です。
　その上、老視になると、英文の活字の小さいこと（笑）、で苦戦しています。
　メガネを二つ掛け替えないとなると、いかに学問がやりにくいか、目だけでも
　若返りたい（笑）。
　

**卵の名無しさん** · 2017/05/06(土) 07:14:59.69

良スレage

**卵の名無しさん** · 2017/05/06(土) 08:54:15.14

>>14
Rのいいところはオープンソースってことだな。
パッケージの関数のソースは
（）なしで入力すると、ソースが表示されるから内部計算の様子がわかる。

例えば
install.packages('Epi')
library(Epi)
twoby2
でtwoby2の関数としてのスクリプトが表示される。

ソースを改造して使えるので便利。
といっても
俺の実力じゃ関数の引数を増やすくらいの改造しかできないけど。

...を追加して引数としてcol='lightblue'としてplot関数に与えてグラフの色を変えるとか簡単にできる。

**卵の名無しさん** · 2017/05/06(土) 09:10:25.15

生データなしでのｔ検定結果を再現するのに自作関数でt検定したりすると勉強になる。

T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
}

# A君の彼女は女子大生、B君の彼女は女子高生。
# Y1女子大生n1=100人とY2女子高生n2=100人の胸囲を測定して
# 前者が平均82 , 標準偏差3
# 後者が平均81 , 標準偏差3

T.test(100,100,82,81,3,3)

> T.test(100,100,82,81,3,3)
[1] 0.0193982

シミュレーションデータを作って検証。

> x=82+scale(rnorm(100))*3
> y=81+scale(rnorm(100))*3
> t.test(x,y,var.equal = TRUE)$p.value
[1] 0.0193982

**卵の名無しさん** · 2017/05/06(土) 13:53:44.66

>>17-18
　ラッキー！
　これ、これが欲しかったんです。データの作り方！！
　ありがとうございます！
　

**卵の名無しさん** · 2017/05/06(土) 13:55:06.57

　
　自分のデータがないとまずRとの距離が遠くて、手でデータを10個ぐらい入れても何の感動もなくて
　とにかく感謝。
　

**卵の名無しさん** · 2017/05/06(土) 15:23:48.05

　
　今日も放送大学「心理統計学」を視聴しておりましたが、今日は、
　2つの標本間の差の分析...全部MCMC法で差を出すので正規分布しか
　出てこない。t検定はおまけで説明。
　
　平均の差の確率分布まで出てくるのがおもしろいのですが、用語と手法に
　慣れないと見ていて疲れます...。
　

**卵の名無しさん** · 2017/05/06(土) 16:39:10.74

当直スレに以前書いたクイズ

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)で統計的に有意差ありという医学ネタにすると面白くないね。

シミュレーションでの計算例は
http://egg.2ch.net/test/read.cgi/hosp/1428282054/756-757

**卵の名無しさん** · 2017/05/06(土) 17:10:55.66

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は胸囲の差
分布をグラフにすると
http://i.imgur.com/huoDBIq.jpg

JDは女子大生のバストが女子高生より5cm以上大きい確率、
JKは女子高生のバストが女子大生より5cm以上大きい確率。

新薬の切り替えたときに臨床的有益性がどれくらいあるかの確率が計算できる。

**卵の名無しさん** · 2017/05/06(土) 18:42:59.31

>>22
　今日の放送大学の講義で、pって何だったんですかって話になって、
　　「元の確率分布から非常に離れた値を取る確率」
　という説明で、
　　「元の確率分布が正しい確率が1-pになるわけじゃないですよ」
　と説明されていた。
　有意差検定とベイズ統計学の差が出ているのかなと思いました。
　

**卵の名無しさん** · 2017/05/06(土) 21:48:23.39

>>24
MotuluskyのIntuitive Biostatisticsに

ｐ値がでたら、帰無仮説は何かを考えよ　という趣旨のことが書いてあった。

リスク比RR＝１　（またはリスク差RD＝０）を帰無仮説にしてｐ値を算出するけれど

RR＝１に帰無仮説を設定する必然性はない、例えばRR＝１．５に設定してｐ値が算出できるという記載が

Rothmanの疫学(introductionの方)にあって

そのグラフｐ値関数を再現できるまでずいぶんと悩んだ（寝当直の暇つぶしになった）。

**卵の名無しさん** · 2017/05/06(土) 23:44:38.35

>>24
その誤解がアニメになっている。

https://www.youtube.com/watch?v=ax0tDcFkPic&;feature=youtu.be

p.value functionを検索していたら

http://sphweb.bumc.bu.edu/otlt/mph-modules/ep/ep713_randomerror/EP713_RandomError7.html

に掲載されていた。

**卵の名無しさん** · 2017/05/07(日) 00:46:48.66

>>26
　かわいいアニメですが、日本語のアニメと同時には聞き取れませんね。今日は無理か、
　ガンダム観たいし（笑）、明日にします。
　後、やっぱり、有意性検定で書かれた論文が多いので、便利なベイズ統計と二足の
　ワラジになりそうですね。
　
　まだ、リスク比とハザード比の差が分からないのですが...。
　

**卵の名無しさん** · 2017/05/07(日) 01:13:01.15

>>27
リスクは一定、
ハザードは時間の経過と共に変化する。

外科治療と内科治療でのハザード比は一定というモデルは
周術期と安定期では一定でないと思う。

リスク比とハザード比をひっくるめてRate Ratioと呼んだりする。
RRはどちらの略か区別がつかないこともしばしば。
比例ハザードモデル=ハザード比一定
ならリスク比と区別がつかなくてもいいんだろうと理解している。

**卵の名無しさん** · 2017/05/07(日) 15:19:34.12

>>28
　つまらない、ピエロみたいな話にレスしてくれてありがとうございます。
　コホート研究と、一回の2x2の検定の差は分かります。

　で、コホート研究の確率分布関数が指数関数という仮定に基づいて
　いませんか？
　それは、統計学者にとっては致命的ではないのでしょうか？
　私は、それががん治療における
　　「最終命題」
　じゃないのかと思うわけです。
　「つらい」の変換で「痛い」がでないんですか？
　「つらい」は日本人的には「痛い」なんですよ。
　
　先生が英語で話したい理由は分かります。世界中に関する影響力が
　違いますよね。
　英語で聞き、英語で表現する、それが確率的に影響力最大なんですよね？

　　「日本人は苦労して、世界に届いている」
　どうしてそれだけの能力があるのに、世界で認められないのか
　そういう話ですよね。
　

**卵の名無しさん** · 2017/05/07(日) 15:22:50.66

　
　確実に臨床統計は、
　　「日本人の生死とQOLを握っている」
　わけでしょ。
　日本人が複数である以上、それを否定する理論があれば、間違いですよ。
　

**卵の名無しさん** · 2017/05/07(日) 15:26:28.75

　
　医師が努力すれば、
　　「医師法第一条に定められた、公衆衛生の利益が守られるわけじゃない」
　という時代に、

　利益最大の医師会や病院団体と、
　厚労省の利益は合わないでしょ？
　

**卵の名無しさん** · 2017/05/08(月) 14:46:04.38

>>30-31
　ごめんなさい、
　他の板の関係で、数学板やら物理板の関係でこんなことを書きましたが、
　
　学問は中立ですよね。中立ゆえの権利があるはずですね。
　と思っているわけです。
　
　数学板でモンティ・ホール問題に疑問を呈して、恣意的な確率操作に対して
　意を表明していますが、
　人間は、理論より納得でしょ。それを得られるかどうか、やってみているわけです。
　ダメなら、呪術をmedicineと呼んでいるアフリカ諸国と同じですよ。
　

**卵の名無しさん** · 2017/05/08(月) 14:52:32.64

　https://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q13120404188
　
　日本語で、「ハザード比、オッズ比、リスク比」を検索すると
　こういうyahoo知恵袋が一番上に出てくるわけで、
　　情けない
　と思うわけです。
　
　誰が日本の臨床統計に対して責任を持っているのか、そういう学会があれば、その検索結果が
　トップにでるように努力すべきなのではないですかね？
　
　ネット時代に対応していない学会は捨てられるか、利益団体に利用されるだけかな？
　

**卵の名無しさん** · 2017/05/09(火) 15:02:02.15

　
　東京馬鹿大で、
　　「医学生が医師を刺した」
　という事象はどう説明するのかね？
　
　明らかに、これは確率論・統計学的な話と言えないかな？
　
　まともな臨床統計を目指すならコメントすべきじゃないし、
　言いたいことがある人なら、コメントすべきじゃないかな。
　
　日本語は使いにくいかな。
　

**卵の名無しさん** · 2017/05/09(火) 15:05:27.95

　
　もっと言えば、
　　「日本人であること、日本国籍を持っていることは有利なのか？」
　って話だと思うよ。
　
　反論したければ、公的保険に対する期待値を示すべきだね。
　

**卵の名無しさん** · 2017/05/13(土) 14:05:24.66

　
　AI=機械学習が人間の決定論に対して、鏡になってくれているわけかな。
　EBMで、医療すら確率論的空間になってしまったが、
　　「人間という機械は、するかしないかの二択でしかあり得ない。」
　だね。
　
　　「確率を提示するのが医師」
　なのかな？　論理学ならすべての事象を集合論で網羅して、真偽を与える
　わけなんだけれど。
　
　私はおかしげな当時の決定論に疑問を持って、大学院は拒否したわけだけれど、
　多くの患者さんの犠牲の下に、私は地位を確立したかったのか、それを拒否したのか、
　その感覚的判断は正しかったのか、バカみたいな論文を書き続けて今の日本の

　医学的な指導者になっているのか、結局、
　　「屍の上に学問は成り立っているわけ」
　でしかない。それをやった奴は手を血に染めているわけだ...。
　その実験をするためには、それを誘導するための意志があったわけだろ。
　
　二重の事後確率を決定するわけだよね。
　

**卵の名無しさん** · 2017/05/13(土) 14:13:03.97

　
　だいたい、普通の臨床医は
　　「即断即決」
　を求められているから。
　
　EBMは嫌いなんじゃないかな、と思う、私もそうだから。
　
　外来の決められた時間で
　　「決定論的結論」
　を求めている患者さんには、この考察にかかる時間は
　　「不利益」
　なんだよね。
　
　確率論は熱力学と同じで、時間的問題を無視している。
　サイコロはいつ投げられても同じ、
　　「コンビニの店長は、時間という価値と確率論的価値
　　はどう判断しているのか？」
　簡単なORにおける在庫管理じゃないよね。
　

**卵の名無しさん** · 2017/05/14(日) 13:05:43.27

>>37
治療による診断はさしずめ、シュレーディンガーの猫かな？

**卵の名無しさん** · 2017/05/14(日) 13:10:25.01

>>37
天気予報にも降水確率つきで発表されるから、正診確率とともに診断を伝えてもいいんじゃないか？

事前確率は往々に経験則もしくは主観的設定だけど。

**卵の名無しさん** · 2017/05/14(日) 16:44:04.43

スレチと思うが、AIは医者業務（手技や承認作業を除く）の何を置き換える事ができるかな？

という問いに対して医者から「絶対これは譲りたくない」と心理的抵抗が強い業務こそ、AIを開発することによって医者減らしに貢献すると思うがどうかな？（救命救急士が行う気管内挿管のように、医者の既得権益？に踏み込む事だろうから。）

**卵の名無しさん** · 2017/05/14(日) 17:56:55.21

>>40
画像診断（ミクロの画像診断、病理も含む）だな。

**卵の名無しさん** · 2017/05/15(月) 00:11:14.10

　
　あの話題の線虫によるがん検診の話がどのぐらいの信憑性があるのか、
　知りたいなぁ。
　

**卵の名無しさん** · 2017/05/15(月) 00:24:43.07

>>40
　AIの能力次第だと思うんだけれど...。

　大量のデータから診断や治療法を探す検索系と、
　機械学習で得られた成果を組み込んだ、画像判断ソフト
　みたいな判断系と、熟練の技術者のデータを入れた
　自動車組み立てロボットみたいな手技系が
　できるんじゃないのかな？
　

**卵の名無しさん** · 2017/05/20(土) 14:35:39.61

　
　土曜の午後は放送大学の無料放送で
　　「ベイズ統計学の復習」
　の予定だったのですが、数週間できませんでした。
　今日は見ます。
　

**卵の名無しさん** · 2017/05/20(土) 14:48:47.08

>>44の続き
　今日も無慈悲に
　　飛ばしている講師の先生。
　不偏分散の話もどこかの回で無視していたのだが、
　　nで割るのとn-1で割るのと
　たいした違いがない。分かるんだけれどね...。

　http://eman-physics.net/math/figures4.html
　にきれいな証明があるからご参考までに。
　

**卵の名無しさん** · 2017/05/20(土) 14:58:43.46

　
　R言語で積み重ねられた、
　　統計学に必要なパッケージと
　　ランダムに作る、100万以上の仮想的なデータ
　があれば、
　統計学の先人達も、そりゃ対数表に頼っていたニュートンやら
　昔の数学・物理学者も1,000円ぐらいの電卓に負けるかもね。
　
　考慮すべき技術的な基盤技術が違う...、
　事前確率分布の関数ことが人間の決められる時代なのかな。
　

**卵の名無しさん** · 2017/05/20(土) 15:13:07.04

　
　あれ、今日もp値の評価の仕方かな。
　最後に有意性検定の否定かな？
　
　まず、帰無仮説の否定...
　まあ、母集団とサンプル群の値が同じだとは言えないよな。
　　サンプルを取る前から偽...
　そりゃあ、恣意的なのかもな。

　Jargon!!!!
　笑います。
　
　ずごい！　刑事事件になってきたぞ。
　アリバイが検定に導入されている。
　　アリバイがないことで有罪の確率が変わるのはおかしい、
　まあそういうことかな。
　

**卵の名無しさん** · 2017/05/20(土) 23:59:15.59

>>47
親子鑑定にベイズが使われたりするんだよな。
理由不十分の原則とかに陪審員が騙されるんだよね。

**卵の名無しさん** · 2017/05/22(月) 18:02:30.07

npoでお金悩み相談。

日々の生活での返済、お支払いでお悩みの方。

急な出費などで、今月の生活費が足りない方。

多重債務、ヤミ金、家賃滞納でお困りの方。

お金に関するお困り事や法的トラブル等HPに記載以外の事でも、お気軽にご相談下さい。

東京、神奈川、千葉、埼玉にお住まいの方は優遇です。

詳しくはHPをご覧下さい。
ｴｽﾃｨｰｴｰで検索

**卵の名無しさん** · 2017/05/24(水) 19:20:21.54

女医が担当医の方が生存率が高いという論文
（１）Comparison of Hospital Mortality and Readmission Rates for Medicare Patients Treated by Male vs Female Physicians
http://scholar.harvard.edu/files/yusuketsugawa/files/tsugawa_jama_im_2016.pdf

若い医師が担当医の方が生存率が高いという論文
（２）Physician age and outcomes in elderly patients in hospital in the US:observational study
http://www.bmj.com/content/357/bmj.j1797

この二つの論文から

（３）若い女医が担当医のときが生存率が最も高い

（４）老いた男性医師が担当医のときが生存率が最も低い

という結論を導きだせるか。

**卵の名無しさん** · 2017/06/10(土) 16:09:52.77

ありゃ、しばらくアクセスしていないうちにスレ落ちちゃったわけで
　申し訳ないです。

確率論におけるヒルベルト空間での議論で
はあ、有意差検定というのは、
　10,000人のデータを真にしたとき
の30例のデータが優位であるかないかの話で、
　10,030例の
検証に至っていないわけです。

臨床統計に至る100例未満の検証が
　10,000の線形空間なのか
　10,100の線形空間なのか
定義がずれているわけで、
AIやらニューロネットワーク的には、後者を選択するわけですよね？
　

**卵の名無しさん** · 2017/06/10(土) 16:25:55.95

　
　明らかにn=10,000がn->無限と規定したいるから
　おかしい話になるわけだ。
　

**卵の名無しさん** · 2017/06/16(金) 15:01:52.10

お久しぶりです、スレ立てしてもメンテができず
　メンテしてくれている人に感謝です。

政治ネタは下げたいのでご容赦を。
　安倍政権のお友達の確率を二項分布でpとしたときに
　明らかに3つ出れば有意な差が出る
でしょうね。

　医学部を増やす、何の意味があるのでしょうか
　竹中さん？
と思いますよ、今の医学系の経営リスクを増やすんでしょうな。

**卵の名無しさん** · 2017/06/16(金) 15:07:56.72

>>50
さて先生のご不満はご承知の上で挑みますよ。
　臨床統計における
　　独立したヒルベルト空間における独立したベクトルをどう設定するかでしょ？
主治医が男性か女性、あるいは患者さんにとって異性か同性なのか、
　それを独立したベクトルとして採用するかどうか
が問題なんでしょうね。
　
臨床統計にとって、
　ヒルベルト空間で何が独立した単位ベクトルなのか
もちろん、集めた母集団で主因子分析をするわけですが
　「しゅいんし」とかな漢字変換に渡したら、「手淫し分析」
が一番に上がるぐらい自己満足感満載です。
　

**卵の名無しさん** · 2017/06/16(金) 15:12:47.34

　
　ここ数日の国会中継を見ていますが、
　　加世学園の問題は確率論的にベイズ更新を行って
　安倍総理の
　　ここ数年の法案提出 → 国会での可決
　確率pがが落ちているのを
　覚悟すべきだと思います。
　
　既にp自体が低下しているのに泥船に乗る確率を
　落とすべきですよね。
　
　ゲーム理論で表現される、
　混合戦略に入りつつあるわけです。
　

**卵の名無しさん** · 2017/06/17(土) 18:17:54.34

中元を配布したリストの提出を求められて税務署に提出。
税務署が「無作為」抽出（実は無作為抽出でなく作為抽出）して調査した５例中、中元を受け取ったのは０であったという。
それをもって税務署は中元は１００例全例虚偽であると認定した。
これはサンプリングに伴うバラつきだと主張して全例への課税を軽減したい。
どういう計算をすればいいか？

**卵の名無しさん** · 2017/06/18(日) 16:18:56.05

# 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
http://i.imgur.com/bhCzo9w.png

**卵の名無しさん** · 2017/06/18(日) 16:40:37.77

当たり確率を連続量の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

http://i.imgur.com/XYCecoG.png

**卵の名無しさん** · 2017/06/19(月) 06:33:53.85

>>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

シミュレーションと少し違うな。

**卵の名無しさん** · 2017/06/19(月) 06:36:19.05

中元を送った期待値はNN/(N+1)でいいかな？

**卵の名無しさん** · 2017/06/22(木) 15:04:26.25

お久さです。はあ、
　世田谷区でデング熱と診断された人が蚊に刺された
として、デング熱の発症確率は増えるんでしょうか？

根拠はあるのですか？　行政的な出費をする根拠になるのですか？

**卵の名無しさん** · 2017/07/01(土) 15:48:51.40

お久しぶりですが、
AIの登場で臨床統計も変わってくると思いますよね。
AIなんてベイズ統計学そのもので検定主義の統計学とは異なるわけです。

死亡率と治療手段を独立とするベクトル空間における内積で議論するには
問題があるように思えますが。

個々の患者さんを説得するには明らかに分散が変わってくると思いますよね？

**卵の名無しさん** · 2017/07/01(土) 20:17:24.70

>>62
ベイズ統計で進めていくときに
説明変数が独立かどうかの設定が難しいと思う。
どうやるのか知らないからではあるが。

模擬試験の判定ってロジスティック回帰の応用だろうけど
英語の点数と国語の点数って独立じゃないだろね。

最近は、ジェネラリストのための内科診断リファレンス
をつまみ読みして生半可な知識を強化している。
明日は雨でしょうから
明日の降水確率は**%
に昇華。

知識が正確になると却って混乱することも増えた。
バセドー病でも抗TPO抗体が陽性になることもあるとか。

**卵の名無しさん** · 2017/07/21(金) 15:51:33.35

ドクターG
http://www4.nhk.or.jp/doctorg/
ありますよね。今でもNHK-Gで毎週放送しています。

新感覚！病名推理エンターテインメント番組「総合診療医ドクターG」。
病名を探り当てるまでの謎解きの面白さをスタジオで展開する！
あなたの症状も解き明かされるか！

というNHK的なクレジットがあるのですが、
　患者さんの生死とかQOKがゲーム的に扱われている
のに疑問を感じます。

診断をどう扱うか、確率のヒルベルト空間では
　ベイズ更新的
な発想で考えていますよね？
一つの事実で1,000-10,000ある診断ラベルのどれの
　確率が変わるのか
という話です。

統計学・AI的思考の話。
　研修医は診断AIに学ぶべき
だと思いませんか？

**卵の名無しさん** · 2017/07/21(金) 16:01:59.58

ある確率密度関数
　P(x)
に関して
　時間tの要素があるんじゃないか
P(x,t)と思いますよね。
　2000年の評価と2017年の評価は違う
としたら、
　2000年に治療を受けていた患者さんと
　2017年に治療を受けていた患者さんは
　別だ
と言いたいわけでしょ？

それで
　サイコロを2回振るのか二つのサイコロを1回振るのか
の違いが分かると思いますよ。

**卵の名無しさん** · 2017/07/21(金) 16:05:13.19

今日、午後4時にNHK-Gにチャンネルを合わせると
　白鵬の土俵入りが観られる
という偶然は偶然じゃなくて
それなりに
　ポテンシャルを変更している
という話が量子論的・確率論的な話だよね。

**卵の名無しさん** · 2017/07/21(金) 16:08:29.61

千秋楽に
　2人の横綱が対決する
ことを事後確率にして
　何人の横綱が必要か
という問題を解いているわけだね。

**卵の名無しさん** · 2017/08/17(木) 20:34:13.43

医者ならば、シリツ卒なら馬鹿である　　から
シリツ卒ならば、医者ならば馬鹿である　が、導けるか？

という論理命題の問題を臨床適用ｗするとこうなる。

甲状腺癌ならば、未分化癌なら予後不良である　　から
未分化癌ならば、甲状腺癌なら予後不良である　が、導けるか？

Rに真理表を使って解答させるスクリプトは以下の通り。

# 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

**卵の名無しさん** · 2017/09/15(金) 02:08:09.73

ここは臨床統計？

**卵の名無しさん** · 2017/09/15(金) 02:58:26.93

>>69
>50できる？

**卵の名無しさん** · 2017/09/15(金) 03:01:19.88

>>70

(3)若い女医の担当患者の成人率が、老いた男医の担当患者の生存率より高いと言えるか？
も追加。

**卵の名無しさん** · 2017/09/24(日) 19:46:12.50

「くじ引きが無作為である」という帰無仮説のもとで宝くじに当選する確率はとても低い（０．０５未満)。
宝くじに当選者がでたということはp<0.05のことが起こったので「くじ引きが無作為」という帰無仮説は棄却される。
正しい議論か？

**卵の名無しさん** · 2017/09/24(日) 20:21:22.18

>>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] #最頻値

**卵の名無しさん** · 2017/09/24(日) 22:09:52.30

> 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

**卵の名無しさん** · 2017/09/24(日) 22:53:55.42

# ド底辺特殊シリツ医大は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))

**卵の名無しさん** · 2017/09/24(日) 23:26:12.20

期待値　sum(pp*f(pp)/sum(f(pp)))

n=600
r=20
n　　　　　　　　　　　　ｎ
Σj*(n-j)Ｃr/nＣr ÷ Σ(n-i)Ｃr/nＣr
j=0　　　　　　　　　　　i=0

**卵の名無しさん** · 2017/09/25(月) 01:09:04.22

## 全体.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)

**卵の名無しさん** · 2017/09/25(月) 06:39:09.92

# 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)

**卵の名無しさん** · 2017/09/25(月) 06:40:55.63

# 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)

**卵の名無しさん** · 2017/09/25(月) 07:58:57.17

##
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)

**卵の名無しさん** · 2017/09/25(月) 09:02:50.99

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)

**卵の名無しさん** · 2017/09/25(月) 15:08:29.67

>>156
素数Aiを任意の数ｎ個集めて
n
Ｂ = Π Ai　 +　 1
i=1

とすると

Ｂを割り切るAi以外の素数が存在することになり、ｎ個集めればn+1個以上の素数の存在が示せることになる。

例

2*3+1で素数7
3*5+1で素数2
3*7+１で素数2と11

**卵の名無しさん** · 2017/09/25(月) 20:13:55.84

# 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)

**卵の名無しさん** · 2017/09/25(月) 20:32:22.48

# 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)

**卵の名無しさん** · 2017/09/25(月) 20:35:29.86

> 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
>

**卵の名無しさん** · 2017/09/25(月) 20:38:01.19

# 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)

**卵の名無しさん** · 2017/09/25(月) 20:48:49.89

>83が可読性が一番いいな

**卵の名無しさん** · 2017/09/26(火) 10:38:08.14

# 全体.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)
}

**卵の名無しさん** · 2017/09/26(火) 10:38:34.73

##
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)

**卵の名無しさん** · 2017/09/26(火) 15:09:39.24

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)

**卵の名無しさん** · 2017/09/26(火) 19:18:43.02

> 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

**卵の名無しさん** · 2017/09/26(火) 20:12:10.07

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))

**卵の名無しさん** · 2017/09/26(火) 20:25:43.72

> 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

**卵の名無しさん** · 2017/09/27(水) 08:43:44.49

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)

**卵の名無しさん** · 2017/09/27(水) 08:44:44.68

> atari.Q(100,10,2)
2.5% 50% 97.5%
6 23 50
mean mode
24.5 20.0
>

**卵の名無しさん** · 2017/09/27(水) 10:05:45.43

> 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

**卵の名無しさん** · 2017/09/27(水) 16:26:42.51

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)

**卵の名無しさん** · 2017/09/27(水) 16:27:31.48

>>97
シリツ医大進学予備校"どていへん予備校"では５年連続不合格、
別のシリツ医大進学予備校"うらぐち予備校"では１０年めで初めて合格者がでたという実績があるとき、
どていへん予備校　と　うらぐち予備校ではどちらが合格可能性が高いと言えるか？

**innuendo** ◆kCkk5BVA12 · 2017/09/27(水) 16:37:18.00

巨乳女子大100人と桃尻女子大100人のどちらが誘いに応じやすいか検討してみる。
巨乳女子大では10人に声をかけて2人が誘いに応じた。
桃尻女子大では無作為順に声をかけていったら5人めが誘いに応じた。
どちらが誘いに応じやすいといえるか？

http://i.imgur.com/DgsfuyX.png

> 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
>

**卵の名無しさん** · 2017/09/27(水) 17:32:04.45

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)

**卵の名無しさん** · 2017/09/27(水) 17:40:20.17

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'))

**innuendo** ◆kCkk5BVA12 · 2017/09/27(水) 17:59:40.63

http://i.imgur.com/DgsfuyX.png

平均値と信頼区間のどちらを選択するかだな。

**卵の名無しさん** · 2017/09/27(水) 20:59:53.74

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)

**卵の名無しさん** · 2017/09/27(水) 21:27:00.63

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)

**卵の名無しさん** · 2017/09/27(水) 23:32:37.60

> 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

**卵の名無しさん** · 2017/09/28(木) 00:10:25.54

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)

**卵の名無しさん** · 2017/09/29(金) 12:14:08.07

1 st 2 nd 3 rd 4 th と表示させるための1行スクリプト

th=ifelse(0<r&&r<4,c('st','nd','rd')[which(r==1:3)],'th')

**卵の名無しさん** · 2017/09/29(金) 19:11:58.46

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')

**卵の名無しさん** · 2017/09/30(土) 18:41:04.08

# 元利均等返済毎月の返済額は一定
.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)

**卵の名無しさん** · 2017/09/30(土) 18:41:31.88

## 元利均等返済 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)
}

**卵の名無しさん** · 2017/10/01(日) 12:38:58.98

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))
}

**卵の名無しさん** · 2017/10/01(日) 12:50:45.81

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)))
}

**卵の名無しさん** · 2017/10/01(日) 12:51:32.13

#
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)

**卵の名無しさん** · 2017/10/01(日) 17:13:29.79

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))

**卵の名無しさん** · 2017/10/03(火) 10:21:03.41

# 百発百中
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)

**卵の名無しさん** · 2017/10/03(火) 14:56:19.98

# ｎ発ｎ中　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)

**卵の名無しさん** · 2017/10/03(火) 14:56:47.23

##
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='全発命中割合')

**卵の名無しさん** · 2017/10/04(水) 11:34:22.68

2-way Contingency Table Analysis http://statpages.info/ctab2x2.html

Epi::twoby2

**卵の名無しさん** · 2017/10/04(水) 13:13:51.18

# 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)
}

**innuendo** ◆kCkk5BVA12 · 2017/10/04(水) 20:50:12.19

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.

**innuendo** ◆kCkk5BVA12 · 2017/10/04(水) 21:03:03.78

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

**innuendo** ◆kCkk5BVA12 · 2017/10/05(木) 11:09:02.87

ゴルゴ13は100発100中
ゴルゴ14は10発10中
ゴルゴ15は1発1中
とする。
各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？

**卵の名無しさん** · 2017/10/05(木) 13:06:15.38

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.

**卵の名無しさん** · 2017/10/05(木) 15:45:07.34

# 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)

**卵の名無しさん** · 2017/10/07(土) 05:54:48.32

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)

**卵の名無しさん** · 2017/10/07(土) 06:05:26.39

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)

**卵の名無しさん** · 2017/10/07(土) 16:41:37.63

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)

**卵の名無しさん** · 2017/10/07(土) 19:50:28.87

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)

**卵の名無しさん** · 2017/10/08(日) 06:36:02.71

>>128
パチンコのスレにこういう記載があった。

>バカは何故かパチンコで勝てると思い込んでいる
>バカは本来なら勝てるのに遠隔や不正で負けていると思い込んでいる
>バカは10分の1で10円が当たるクジより1万分の1で8000円が当たるクジの方が儲かると思ってる
>算数すらできないバカが必死に守って支えてきたのがパチンコ業界

これを読んでこんな問題を考えてみた。

宝くじAは1000本中100本が当たりで当たりは 100万円の賞金、

宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、

どちらも売り出し価格は同じなので100本買うことにする。

どちらを何本買うときに賞金の期待値が最大になるか、シミュレーションしてみよ。

**卵の名無しさん** · 2017/10/08(日) 12:43:27.27

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')

**卵の名無しさん** · 2017/10/08(日) 18:38:07.36

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')

**卵の名無しさん** · 2017/10/09(月) 08:02:11.41

# 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)

**卵の名無しさん** · 2017/10/09(月) 11:07:15.42

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
と動作してくれた。

**卵の名無しさん** · 2017/10/09(月) 14:34:35.08

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)

**innuendo** ◆kCkk5BVA12 · 2017/10/09(月) 15:13:49.58

>>127
http://imagizer.imageshack.com/img923/1874/4NP9BA.png

**卵の名無しさん** · 2017/10/09(月) 15:39:38.20

>>135
http://i.imgur.com/094oClD.png

**卵の名無しさん** · 2017/10/09(月) 15:40:36.46

> summary(g13)
Min. 1st Qu. Median Mean 3rd Qu. Max.
8992 9866 9933 9903 9972 10000
> which.max(table(g13))
10000
714

**卵の名無しさん** · 2017/10/09(月) 19:26:11.72

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(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))

**卵の名無しさん** · 2017/10/10(火) 06:51:19.25

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()

**卵の名無しさん** · 2017/10/10(火) 16:36:50.99

> '+'(1,2)
[1] 3

> '*'(2,3)
[1] 6

> a=1:10

> '['(a,7)
[1] 7

'[' って関数なんだな。

**卵の名無しさん** · 2017/10/10(火) 21:23:23.64

>>129

宝くじAは1000本中300本が当たりで当たりは 100万円の賞金、

宝くじBは1000本中 10本が当たりで当たりは1000万円の賞金、

とAの期待値＝Bの期待値の3倍のとき

http://i.imgur.com/6CN1yOt.png

**卵の名無しさん** · 2017/10/11(水) 06:48:41.73

>>139
続きのコードが書き込めないので画像でアップ
http://imagizer.imageshack.com/img922/5007/LKoY9v.jpg

**卵の名無しさん** · 2017/10/11(水) 07:14:30.46

# 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)

**卵の名無しさん** · 2017/10/11(水) 20:53:31.06

確かにマジックだな

##-- 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)

**卵の名無しさん** · 2017/10/11(水) 21:31:48.49

[　って関数だったんだ

(d=outer(11:19,11:19,'*'))
d[,3]
apply(d,1,'[',3)

**卵の名無しさん** · 2017/10/11(水) 23:59:07.55

(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))

**卵の名無しさん** · 2017/10/12(木) 18:41:56.42

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

**卵の名無しさん** · 2017/10/14(土) 04:34:05.51

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,'*')

**卵の名無しさん** · 2017/10/14(土) 06:39:17.90

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())

**卵の名無しさん** · 2017/10/14(土) 19:44:47.68

##
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))

**卵の名無しさん** · 2017/10/14(土) 19:45:01.23

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)

**卵の名無しさん** · 2017/10/14(土) 19:45:34.01

## 和で直積
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])

**卵の名無しさん** · 2017/10/14(土) 19:46:23.30

> 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

同型写像確認

**卵の名無しさん** · 2017/10/14(土) 22:27:16.20

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)

**卵の名無しさん** · 2017/10/15(日) 07:15:02.95

https://en.wikipedia.org/wiki/Euler%27s_totient_function

を参考にオイラー関数をグラフにしてみた。

http://i.imgur.com/9db3JHt.png

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)

**卵の名無しさん** · 2017/10/16(月) 14:46:14.17

(p+1)^(p^(n-1)) ≡1 (mod p^n)

(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))

**卵の名無しさん** · 2017/10/16(月) 18:24:23.14

# (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)

**卵の名無しさん** · 2017/10/16(月) 19:48:05.23

(p+1)^(p^(n-1)) ≡ p^n + 1 (mod p^(n+1))
p:奇素数

**卵の名無しさん** · 2017/10/18(水) 20:51:15.46

# (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)
}

**卵の名無しさん** · 2017/10/19(木) 15:04:12.35

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]
}

**卵の名無しさん** · 2017/10/19(木) 15:04:48.68

.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

**卵の名無しさん** · 2017/10/19(木) 16:00:22.26

結合側の検証

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

**卵の名無しさん** · 2017/10/19(木) 19:20:37.10

>>162
結合則

**卵の名無しさん** · 2017/10/19(木) 19:26:48.29

演算が閉じている、単位元が存在する、逆元が存在する、の確認は容易だけど

結合則の確認は手計算だと手間がかかる。

プログラムを組んで確認できて納得。

**卵の名無しさん** · 2017/10/20(金) 13:57:10.48

紛らわしいなぁ

群論における「位数」とは、二つの異なる定義があり、それぞれ、
①有限群Gの位数＝有限群Gの元の個数、
②有限群Gの元aの位数＝有限群Gの元aにおいて、a^m=e となる最小の正の整数、
となります。

**卵の名無しさん** · 2017/10/20(金) 14:07:59.55

## 位数２ｎの二面体群の演算表を作る
# 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

**卵の名無しさん** · 2017/10/20(金) 14:08:06.73

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)

**卵の名無しさん** · 2017/10/20(金) 14:08:34.26

> 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

**卵の名無しさん** · 2017/10/20(金) 14:21:16.31

これが抜けていた。

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]
}

**卵の名無しさん** · 2017/10/20(金) 14:30:53.07

## 位数２ｎの二面体群の演算表を作る
# 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

**卵の名無しさん** · 2017/10/20(金) 14:32:25.05

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

**卵の名無しさん** · 2017/10/20(金) 14:32:49.69

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)

**卵の名無しさん** · 2017/10/20(金) 14:41:11.07

>170-172で動作確認

**卵の名無しさん** · 2017/10/21(土) 07:30:34.81

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

**卵の名無しさん** · 2017/10/21(土) 19:00:23.36

## 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)

**卵の名無しさん** · 2017/10/21(土) 19:00:49.55

> (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

**卵の名無しさん** · 2017/10/21(土) 20:08:08.06

## 二面体群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]
}
}

**卵の名無しさん** · 2017/10/21(土) 20:19:32.79

http://i.imgur.com/fgx1jxd.png

**卵の名無しさん** · 2017/10/21(土) 20:20:03.91

##
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)

**卵の名無しさん** · 2017/10/21(土) 20:21:10.03

> (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

# 御明算

**卵の名無しさん** · 2017/10/21(土) 21:11:33.21

>>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)

**卵の名無しさん** · 2017/10/21(土) 23:55:34.83

## 二面体群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

**卵の名無しさん** · 2017/10/21(土) 23:56:01.41

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)
}
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

**卵の名無しさん** · 2017/10/21(土) 23:56:34.66

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

**卵の名無しさん** · 2017/10/21(土) 23:57:03.96

.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)
}

**卵の名無しさん** · 2017/10/21(土) 23:57:45.77

> 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

2017/10/22(日) 09:18:14.91

# 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)

2017/10/22(日) 13:48:21.69

> 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

2017/10/22(日) 15:51:11.54

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

2017/10/22(日) 16:46:17.03

> 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

2017/10/22(日) 17:29:41.55

# 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 ...
フェルマーの定理

2017/10/22(日) 18:21:14.08

>>191
フェルマーも小定理の方ね

**卵の名無しさん** · 2017/10/25(水) 19:11:29.66

>>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]

**卵の名無しさん** · 2017/10/25(水) 19:13:22.53

>>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

**卵の名無しさん** · 2017/10/25(水) 19:13:48.14

> 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

**卵の名無しさん** · 2017/10/25(水) 19:24:32.39

>>195
n=10　の時のPDFと９５％CI下限値0.7615958

https://i.imgur.com/VrACoar.png

**卵の名無しさん** · 2017/10/25(水) 19:24:57.43

>>195
n=10　の時のPDFと９５％CI下限値0.7615958

http://i.imgur.com/VrACoar.png

**卵の名無しさん** · 2017/10/25(水) 20:40:32.01

#.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)
}

**卵の名無しさん** · 2017/10/26(木) 12:27:43.35

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)
}

**卵の名無しさん** · 2017/10/26(木) 16:21:26.13

# ｆ回失敗した後にｓ回目の成功,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)
}