0176132人目の素数さん
2018/08/14(火) 22:15:51.04ID:1M/cP5sU省エネ化してRに移植
# 確率を行列化
rm(list=ls())
X = function(n,r=10,w=90){
# s[i,j]
rw=r+w # 試行前総玉数
J=rw+n # jの上限
s=matrix(0,nrow=n,ncol=J) # i回試行後に総数がj個である確率の行列
s[1,rw-1]=w/rw ; s[1,rw+1]=r/rw # 1回試行後
if(n > 1){
for(i in 2:n){
for(j in r:J){ # jはr未満にはならない
# if(j==1) s[i,j] = s[i-1,j+1]*(j+1-r)/(j+1)
if(j==J) s[i,j] = s[i-1,j-1] * r/(j-1)
else s[i,j] = s[i-1,j-1] * r/(j-1) + s[i-1,j+1]*(j+1-r)/(j+1)
}
}
}
total=sum((r:J)*s[n,r:J])
white=total-r
return(c(total=total,white=white))
}
X(2)
546621/5555
vX=Vectorize(X)
vX(1000:1010)
> vX=Vectorize(X)
> vX(1000:1010)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
total 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5 20.5
white 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5 10.5