蒙特卡洛方法求圆面积的R代码

snow

有人问了一个蒙特卡洛求解问题，为了帮助理解蒙特卡洛方法思想，做了个例子，希望有所帮助。

问题：半径是1的圆，写个程序，不许用pi，结果接近3.14

思想就是：
半径为1的圆可以放到边长为2的正方形里面
边长为2的正方形面积是（2*圆的半径)^2
我随机在正方形里面均匀选10000个点，看看有多少个点可以落在圆里，也就是圆面积和正方形面积的比例就算出来了


如果手动试验可以这样做：找张边长为2的纸，然后上面画个内切圆，然后拿一粒沙子，往纸上扔，扔1000次，查查多少次落到了圆内，就把pi求出来了

求解代码如下：

1. area<-function(radius=1,MonteCarlo_point_num=1000){
   
2. randpoint_x=runif(MonteCarlo_point_num,0,2*radius)
   
3. randpoint_y=runif(MonteCarlo_point_num,0,2*radius)
   
4. point_num_in_circle=0
   
5. for (i in 1:MonteCarlo_point_num){
   
6. distance_randpoint_circlecenter=sqrt((randpoint_x[i]-radius)^2+(randpoint_y[i]-radius)^2)
   
7. if(distance_randpoint_circlecenter<radius){
   
8. point_num_in_circle=point_num_in_circle+1
   
9. }
   
10. }
    
11. percent=point_num_in_circle/MonteCarlo_point_num
    
12. area=percent*(2*radius)^2
    
13. return(area)
    
14. }
    
15. area(1,100000)

复制代码

biostat

改进下运算速度

1. area<-function(radius=1,MonteCarlo_point_num=1000){
   
2. randpoint_x=runif(MonteCarlo_point_num,0,2*radius)
   
3. randpoint_y=runif(MonteCarlo_point_num,0,2*radius)
   
4. distance_randpoint_circlecenter<-sapply(1:MonteCarlo_point_num,function(x) sqrt((randpoint_x[x]-radius)^2+(randpoint_y[x]-radius)^2))
   
5. point_num_in_circle<-length(distance_randpoint_circlecenter[distance_randpoint_circlecenter<radius])
   
6. percent=point_num_in_circle/MonteCarlo_point_num
   
7. area=percent*(2*radius)^2
   
8. return(area)
   
9. }
   
10. area(1,100000)
    

复制代码