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
214:卵の名無しさん
17/10/30 07:00:24.65 OjMfeTM6.net
## 写像を元とする集合の積を計算する
equal <- function(x,y) sum((x-y)^2)==0 & length(x)==length(y)
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)
}
witch <- function(x,L=G){
re=NULL
for(i in 1:length(L)){
if(equal(L[[i]],x)) return(i) # i for names[i]
}
}
tikan3 <- function(X,Y,L=G){ # replace Y first then X and return index of L
Z=tikan2(X,Y)
witch(Z,L)
}
product <- function(A,B,L=G){ # product of sets, retur index of L
la=length(A)
lb=length(B
215:) gr=expand.grid(1:la,1:lb) f<-function(x,y) tikan3(A[[x]],B[[y]],L) re=mapply(f,gr[,1],gr[,2]) return(sort(unique(re))) }
216:卵の名無しさん
17/10/30 07:01:59.31 OjMfeTM6.net
GROUP <- function(n=6){
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
return(G)
}
217:卵の名無しさん
17/10/31 16:10:43.91 1m4iMpv8.net
HDCI <- function(PMF,cl=0.95){ # Highest Density Confidence Interval
PDF=PMF/sum(PMF)
rsPDF=rev(sort(PDF))
min.density=rsPDF[min(which(cumsum(rsPDF)>=cl))]
index=which(PDF>=min.density)
data.frame(lower.idx=round(min(index)),upper.idx=round(max(index)),actual.CI=sum(PDF[index]))
}
218:卵の名無しさん
17/10/31 16:11:33.50 1m4iMpv8.net
# n発r中の期待値
Golgo2 <- function(n=3,r=1,cl=0.95,k=0.00001,Print=FALSE){
PMF <- function(x) choose(n,r)*x^r*(1-x)^(n-r)
xx=seq(0,1,by=k)
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=5)
if(Print){plot(xx,pmf,lwd=1,xlab='probability of hit',
ylab=paste('probabilty of',r,'hits out of ',n,'shoots'))}
}
219:卵の名無しさん
17/10/31 20:38:26.69 A7oOP/mA.net
# １年：進級失敗１０人、うち１人放校
# ２年：進級失敗１６人、放校なし
# ３年：進級失敗３４人、うち放校９人
# ４年：進級失敗９人、うち放校２人
# ５年：進級失敗１０人、うち放校１人
# ６年：卒業失敗２６人、うち放校１人
# 一学年約120～130人前後。
#ｽﾚﾘﾝｸ(doctor板:1番)
flunk=c(10,16,34,9,10,26)
expel=c(1,0,9,2,1,1)
n=125
total.fl=sum(flunk)
total.ex=sum(expel)
re=matrix(NA,6,1)
for(i in 1:6){
re[i,]=poisson.test(c(flunk[i],n) ,c(total.fl,n*6))$p.value
}
re
for(i in 1:6){
re[i,]=poisson.test(c(expel[i],n) ,c(total.ex,n*6))$p.value
}
re
prop.test(flunk,rep(n,6))
pairwise.prop.test(flunk,rep(n,6),'bon')
prop.test(expel,rep(n,6))
d=cbind(expel,rep(n,6)-expel) ; d
fisher.test(d)
pairwise.prop.test(expel,rep(n,6),'none')
pairwise.prop.test(expel,rep(n,6),'bon')
220:卵の名無しさん
17/11/01 09:28:50.55 zVXhNw2e.net
## 確率分布から信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=TRUE){
xx=seq(0,1,by=k)
xx=xx[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
print(c(lower=lower,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
HDCI2(function(x)dnbinom(5-1,1,x))
221:卵の名無しさん
17/11/01 09:38:36.29 zVXhNw2e.net
## 確率分布から最尤値・期待値・信頼区間を出す
HDCI2 <- function(PMF,cl=0.95,k=0.0001,Print=FALSE){
pp=seq(0,1,by=k)
xx=pp[-1]
pmf=sapply(xx,PMF)
pdf=pmf/sum(pmf)
rspdf=rev(sort(pdf))
min.density=rspdf[min(which(cumsum(rspdf)>=cl))]
index=which(pdf>=min.density)
lower.idx=round(min(index))
upper.idx=round(max(index))
lower=xx[lower.idx]
upper=xx[upper.idx]
actual.CI=sum(pdf[index])
mpdf=sum(xx*pdf)
mode=xx[which.max(pmf)]
print(c(lower=lower,mode=mode,mean=mpdf,upper=upper,actual.CI=actual.CI),digits=4)
if(Print) plot(xx,pmf,xlab='prior',ylab='posterior')
}
222:卵の名無しさん
17/11/05 17:08:17.51 G4ZpDCNG.net
NHST Has 100% False Alarm Rate in Sequential
Testing
Under NHST, sequential testing of data generated from the null
hypothesis will eventually lead to a false alarm. With infinite
patience, there is 100% probability of falsely rejecting the null.
This is known as “sampling to reach a foregone conclusion” (e.g.,
Anscombe, 1954). To illustrate this phenomenon, a computer
simulation generated random values from a normal distribution
with mean zero and standard deviation one, assigning each sequential
value alternately to one or the other of two groups, and at each
step conducting a two-group t test assuming the current sample
sizes were fixed in advance. Each simulat
223:ed sequence began with N1 N2 3. If at any step the t test indicated p .05, the sequence was stopped and the total N ( N1 N2) was recorded. Otherwise, another random value was sampled from the zero-mean normal and included in the smaller group, and a new t test was conducted. For purposes of illustration, each sequence was limited to a maximum sample size of N 5,000. The simulation ran 1,000 sequences. NHST <- function(N=5000){ x=rnorm(3) y=rnorm(3) while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){ x=append(x,rnorm(1)) y=append(y,rnorm(1)) } return(length(x)) } res <- replicate(1000,NHST())
224:卵の名無しさん
17/11/05 17:27:50.94 G4ZpDCNG.net
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
# Under NHST, sequential testing of data generated from the null
# hypothesis will eventually lead to a false alarm. With infinite
# patience, there is 100% probability of falsely rejecting the null.
# This is known as “sampling to reach a foregone conclusion” (e.g.,
# Anscombe, 1954). To illustrate this phenomenon, a computer
# a computer simulation generated random values from a normal distribution
# with mean zero and standard deviation one, assigning each sequential
# value alternately to one or the other of two groups, and at each
# step conducting a two-group t test assuming the current sample
# sizes were fixed in advance. Each simulated sequence began with
# N1=N2=3. If at any step the t test indicated p .05, the
# sequence was stopped and the total N (=N1+N2) was recorded.
NHST <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & t.test(x,y,var.equal=TRUE)$p.value >= 0.05){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120,10240)
RjR=sapply(NN,FRR)
plot(NN,RjR,type='b')
225:卵の名無しさん
17/11/05 18:24:17.31 G4ZpDCNG.net
# Bayesian decision making, using the HDI and ROPE, does not
# suffer a 100% false alarm rate in sequential testing. Instead, the
# false alarm rate asymptotes at a much lower level, depending on
# the choice of ROPE. For illustration, again a computer simulation
# generated random values from a normal distribution with mean of
# zero and standard deviation of one, assigning each sequential value
# alternately to one or the other of two groups but at each step
# conducting a Bayesian analysis and checking whether the 95%
# HDI completely excluded or was contained within a ROPE from 0.15 to 0.15
# ROPE* region of practical equivalence
librry(BEST)
BA <- function(N){
mc=BEST::BESTmcmc(rnorm(N),rnorm(N),numSavedSteps=1000, burnInSteps=500)
re=summary(mc,ROPEm=c(-0.15,15))
return(re[['muDiff','%InROPE']])
}
NN=c(10,20,40,80,160,320,640,1280,2560,5120)
inrope=sapply(NN,BA)
plot(NN,inrope,type='b',pch=19,xlab='N',ylab='% in ROPE')
226:卵の名無しさん
17/11/05 20:18:58.81 G4ZpDCNG.net
Bayesian Estimation Supersedes the t-test (BEST) - online
URLﾘﾝｸ(www.sumsar.net)
トレースプロットがみられるのがうれしい。
227:卵の名無しさん
17/11/05 20:22:56.18 G4ZpDCNG.net
新版の予約受付中か。
URLﾘﾝｸ(www.amazon.co.jp)
ベイズ統計がどれくらい組み込まれたをみてから買うかな。
228:卵の名無しさん
17/11/05 20:28:14.21 G4ZpDCNG.net
# A君の彼女は女子大生、B君の彼女は女子高生。
# Y1女子大生n1=100人とY2女子高生n2=100人の胸囲を測定して
# 前者が平均82 , 標準偏差3
# 後者が平均81 , 標準偏差3
# 有意差はあるか？
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)
}
T.test(100,100,82,81,3,3)
y1=82+scale(rnorm(100))*3
y2=81+scale(rnorm(100))*3
t.test(y1,y2,var.equal = TRUE)
library(BEST)
BESTout <- BESTmcmc(y1,y2)
plot(BESTout,ROPE=c(-2,2))
summary(BESTout,ROPEm=c(-2,2))
plotAll(BESTout,ROPEm=c(-2,2))
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-0.5,0.5), ROPEsd=c(-1,1), ROPEeff=c(-0.5,0.5),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=5)
powerPro
229:卵の名無しさん
17/11/05 20:33:01.04 G4ZpDCNG.net
ProData <- makeData(mu1=82, sd1=3, mu2=81, sd2=3, nPerGrp=100,
pcntOut=10, sdOutMult=2.0, rnd.seed=NULL,showPlot=TRUE)
proMCMC <- BESTmcmc(proData$y1, proData$y2, numSavedSteps=2000)
N1plan <- N2plan <- 100
powerPro <- BESTpower(proMCMC, N1=N1plan, N2=N2plan,
ROPEm=c(-2,2), ROPEsd=c(-0,0), ROPEeff=c(-0,0),
maxHDIWm=5.0, maxHDIWsd=2.0, maxHDIWeff=1.0, nRep=1000)
powerPro
230:卵の名無しさん
17/11/05 20:51:07.47 G4ZpDCNG.net
Jikkan <- function(x){
re=summary(BESTout,ROPEm=c(-x,x))
re[['muDiff','%InROPE']]
}
xx=seq(0,2.5,by=0.01)
InRope=sapply(xx,Jikkan)
plot(xx,InRope,type='l',lwd=2, xlab='Difference',ylab='% Practically Equivalent')
abline(h=95, lty=3)
(Diff=uniroot(function(x,u0=95) Jikkan(x)-u0, c(1.5,2))$root)
plot(BESTout,ROPE=c(-Diff,Diff))
URLﾘﾝｸ(i.imgur.com)
231:卵の名無しさん
17/11/06 21:07:39.66 1uZMaFLW.net
##
library(BayesFactor)
tBF <- function(x,y){
t=t.test(x,y,var.equal = TRUE)$statistic
ttest.tstat(t,length(x),length(y),simple=TRUE)
}
NHST2 <- function(N){
x=rnorm(3)
y=rnorm(3)
while(length(x) < N & tBF(x,y) <= 1){
x=append(x,rnorm(1))
y=append(y,rnorm(1))
}
return(length(x))
}
FRR2 <- function(N,k=100){ #False Reject Rate
re <- replicate(k,NHST2(N))
return(mean(re!=N)) # rate falsely rejected as significantly different
}
> FRR2(10)
[1] 0.37
> FRR2(50)
[1] 0.46
> FRR2(100)
[1] 0.56
ベイズ因子をつかっても ｔ 検定を組み込むと
# NHST Has 100% False Alarm Rate in Sequential Testing (NHST : Null Hypothesis Significance Testing)
という結論は変わらないな。