20/07/14 16:54:04.59 bf2JlY+k.net
# 分数表示
library(gmp)
pmf <- function(n.taxi,Max=60){
p=list() # 確率の配列
for(n.obs in 1:Max){ # n.obs : 観察台数1 ~ 60
# 総台数n.taxiで観察数がn.obsのとき最大番号がMaxになる確率
p[[n.obs]]=chooseZ(Max-1,n.obs-1)/chooseZ(n.taxi,n.obs)
}
s=as.bigq(0)
for(i in 1:Max) s = s + p[[i]] # sに分数加算
s/Max # 平均値(期待値)
}
n=60:100 # 総台数の候補
pdf=pmf(n)/sum.bigq(pmf(n)) # pmfをpdf化
(E=sum(n*pdf)) # 期待値
asNumeric(E)
# 観察された台数が不明なときのシミュレーション
sim <- function(){
M=m=0 # m:観察台数 M:最大番号 (初期値0)
while(M!=60){ # M=60でないなら
N=sample(60:100,1) # タクシー総数Nを60 ~ 100から選ぶ
m=sample(1:min(N,60),1)# 観察する台数mを1 ~ min(N,60)から選ぶ
M=max(sample(1:N,m)) # N台からm台選択して最大値をMにいれる
}
return(N) # タクシー総数を返す
}
re=replicate(1e4,sim()) # 1万回繰り返して平均値(期待値)を算出
mean(re)
BEST::plotPost(re,xlab='N',cex.lab=1)