#碰到钉子后的反应函数
fun0301=function(t0,p=0.5){
c(t0[1]+1,t0[2]+ if(runif(1)>=p) 1 else -1 )
}
#画钉板的函数
fun.up <- function (n) {
plot(0,type="n",bty="n",xaxt="n",yaxt="n",
xlim=c(0,2*n+1),ylim=c(1,n+1)
)
axis(2,1:n,n:1)
for(i in 0:(2*n+1)) lines(c(i,i) ,c(0,n),lty=3 )
for(i in0:n ) lines(c(0,2*n+1),c(i,i),lty=3)
points( cbind(rep(c(2*(0:n),2*(0:n)+1),n/2)[1:(n*(n+1))],
rep(n:1,each=n+1) ),
pch=17)
}
#根据钉的位置添加点
fun0303=function(t0,n,...){ points(t0[2],n+1-t0[1]+0.4,...) }
#添加小线段
fun0304=function(t0,t1,n,...){
if(t0[1]<2)stop("所在层数要大于2！")
t2=t0[1]-0:1
lines(t1[t2],n+1-t2+0.4,...)
}
#...
fun03.moive <- function (m=100, n=20, p=0.5,ymax=m/3,sleep02=0.1) {
y=c() #记录球最终的位置
windows(15,15)
layout(matrix(1:2,2) )
for(mi in1:m){
####上面的钉板
par(mar=c(0,2,1,1))
fun.up(n=n)
#球的走动
t0=c(1,n) #第1层，第n个钉的上方
t1=c(n) #记录各层球的位置
fun0303(t0,n,pch=16,col=2) #显示钉的位置
#走动
while(t0[1]<=n){
t0=fun0301(t0,p) #跳至下一层
t1=c(t1,t0[2]) #保存钉的位置
fun0303(t0,n,pch=16,col=2) #显示钉的位置
fun0304(t0,t1,n,lty=1,col=2,lwd=2) #加路线
}
y=c(y,t1[length(t1)]) #保存球最终的位置
####下面的直方图
par(mar=c(1,2,0,1))
hist(y+0.5,breaks=(-1):(2*n+1)+0.5,xlim=c(0,2*n+1),ylim=c(0,ymax),
xaxt="n",main="",xlab="",ylab="",col=6)
for(i in0:(2*n+1)) lines(c(i,i),c(0,max(table(y))),lty=3) #辅助线
z=table(y)
lines(rep(y[length(y)],2),c(ymax,z[names(z)==y[length(y)]] ),
lty=2,col=2,lwd=2 )
text(1,ymax,paste("p=",p,sep=""),col=4 )
Sys.sleep(sleep02)
}
}
fun03.moive(p=0.5)