R语言-蒲丰投针问题(向量化)

蒲丰投针问题的具体过程及推导可参考百科，此处不再赘述，这里仅给出模拟步骤和实现代码。

蒲丰投针试验的模拟步骤

• 设定试验总次数 $N$N，成功次数$m = 0$m=0，针的长度为$h$h，平行线间的距离为$a$a ($h\leqslant a$h⩽a)
• 产生随机数$X\sim U\left(0,a/2\right)$X∼U(0,a/2)，$\theta\sim U\left(0,\pi\right)$θ∼U(0,π)
• 若$X\leqslant \frac{h}{2}sin\theta$X⩽2h​sinθ，$m=m+1$m=m+1，转上一步
• 若达到试验总次数N，停止实验，计算$\pi$π的估计值$\widehat{\pi }=\frac{2hN}{am}$π=am2hN​

代码实现

• 常规思想

buffon<-function(N,a,h){
  m = 0
  for (i in 1:N){
    x<-runif(1,0,a/2)
    theta<-runif(1,0,pi) %X和theta的产生
    if(x<=h/2*sin(theta)){ %判断条件
      m=m+1
    }
  }
  pin<-2*h*N/(a*m) %pi的估计值
  result<-c('pi估计值'=pin)
  return(result)
}

buffon(10000,2,0.8) %调用

• 向量化思想

buffon<-function(N,a,h){
  m = 0
  if(h<=a){
    x<-runif(N,0,a/2)
    theta<-runif(N,0,pi) %直接产生N个x和theta
    m = sum(x<=h/2*sin(theta)) %满足判别条件时值为1，并求和
    pin<-2*h*N/(a*m)
  }else{
    stop("Error") %h小于等于a的判别
  }
  return(pin)
}

buffon(10000,2,0.8)

避免显示循环，遵循向量化原则，才能高效的使用R语言。