>>140
省エネ化して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