编辑代码

#Buffon's needle problem Using computer simulation
#（1）取一张白纸，在上面画上许多条间距为a的平行线。
#（2） 取一根长度为l（l≤a） 的针，随机地向画有平行直线的
#纸上掷n次，观察针与直线相交的次数，记为m。
#（3）计算针与直线相交的概率：
#【Buffon本人证明了】这个概率是：p=2l/πa (其中π为圆周率）
#接下来设计实验，取l=0.5 a=1，那么π=总投针次数/相交次数
#
# 绘制空白图形
plot(c(0,2),c(0,2),type='n',main='布丰投针实验',xlab='X',ylab='Y')
# 增加平行线
abline(h=0.5)
abline(h=1.5,col='red')
finished <- FALSE
# 【trial为实验次数，cross为交叉次数】
trial <- 0
cross <- 0
while (!finished) {
    # Dist为针的中心距离红线的垂直距离
    # Theta为针的角度
    Dist <- runif(1,min=0,max=1/2)
    Theta <- runif(1,0,pi)
    # central.x为针中心点的横坐标
    # central.y为针中心点的纵坐标
    central.x <- runif(1,0.5,1.5)
    central.y <- Dist +1
    # 计算针两端的坐标
    y1 <- sin(Theta)/4 + central.y
    x1 <- cos(Theta)/4 + central.x
    y2 <- sin(Theta+pi)/4 + central.y
    x2 <- cos(Theta+pi)/4 + central.x
    trial <- trial +1
    # 计数交叉次数
    cross <- cross + ifelse(0.25*sin(Theta)>=Dist,1,0)
    # 绘制针的线型和中心点
    lines(c(x1,x2),c(y1,y2),lty=2)
    points(central.x,central.y,pch=16,col='grey')
    cat('trial=',trial,'cross=',cross,'PI=',trial/cross,'\n')
    #continue?
    input <- readline('stop?')
    # 若输入y，则结束实验
    if (input =='y') finished <- TRUE
}