- 昨日の続き
- こちらやこちらの記事を参考にする
- N個の一様乱数でa未満の乱数がk-1個である確率
- それの微分がN個の一様乱数で第k番目に小さい乱数の値がaである確率
- これをＲで、行列やapply()を使ってコンパクトに計算してプロットする
# カットオフp値はベクトルで与えよう
as<-seq(from=0,to=1,by=0.05)
N<-10 # マーカー数
# N個の一様乱数でa未満の乱数がk-1個である確率
aa<-exp(t(log(as)%*%t(0:N)+log(1-as)%*%t(N:0))+(lgamma(N+1)-lgamma((0:N)+1)-lgamma((N:0)+1)))
# 累積確率にしてプロット
aa2<-1-apply(aa,2,cumsum)
matplot(t(aa2),type="l")
# それの微分がN個の一様乱数で第k番目に小さい乱数の値がaである確率
bb<-t(aa*(-(0:N)))/as + t(aa*(N:0))/(1-as)
# その累積
cc<-apply(bb,1,cumsum)
matplot(t(cc),type="l")
- これは、aが小さいときには、アーラン分布とか、ガンマ分布とか(関連項目)？？？要、確認。
