【Buffon蒲丰投针实验】R语言代码：实现投针实验模型可视化

萌新一个，概率论作业，模拟Buffon蒲丰投针实验~~

网上找了一圈，都找不到用R语言做这个实验的代码和教程，于是花了前后一两个月时间，前三天花了三天三夜完成代码编写，参考各种资料，自己编写代码，终于调试成功啦！！后面的时间就是逐步完善代码。已经优化了6版了。

R语言代码（原始版本）展示（去掉函数之后的代码）

这里使用Rstudio来编写代码

注意：这里因为是原始版本，没有使用到函数调参，以及与函数版本运用的判断线段与直线相交的条件不同，所以这里的代码只考虑了投针次数的变化，所以：平行线间距、针的长度、平行线的数目，这三个的值不可更改，只可更改投针的次数times

LineSpacing=1#平行线间距
Needle_length=LineSpacing/2#针的长度
begin_time<-1#第一次投针
times<-10000#投针的次数
alltimes<-begin_time:times#表示从1到投针的次数的向量
count_NeedleIntersect=0#初始化针与线的相交次数

#以下为创建向量,即预分配内存
NeedleX_1=c()
NeedleY_1=c()
NeedleLine_angle=c()
NeedleX_2=c()
NeedleY_2=c()
my_pi=c()

for (i in alltimes)
{
  NeedleX_1[i]=8*runif(1) #随机生成投针的起始点横坐标（乘上8为扩大随机0到1数为0到8）
  NeedleY_1[i]=8*runif(1)#随机生成投针的起始点纵坐标（乘上8为扩大随机0到1数为0到8）
  NeedleLine_angle[i]=pi*runif(1)#随机生成针与线的夹角
  NeedleX_2[i]=NeedleX_1[i]+cos(NeedleLine_angle[i])*Needle_length#计算投针的另一端点横坐标（三角函数）
  NeedleY_2[i]=NeedleY_1[i]+sin(NeedleLine_angle[i])*Needle_length#计算投针的另一端点纵坐标（三角函数）
}
NeedleX<-c(NeedleX_1,NeedleX_2)#合并起始点横坐标和终点横坐标
NeedleY<-c(NeedleY_1,NeedleY_2)#合并起始点纵坐标和终点纵坐标
par(mfrow=c(1,2))##将绘图窗口划为1行两列，同时显示两图
plot(NeedleX,NeedleY,pch=20)#绘制点图
abline(h=0:8,col="red")#绘制平行线
for (i in alltimes)
{
  segments(NeedleX_1[i],NeedleY_1[i], NeedleX_2[i], NeedleY_2[i],col = "blue")#连接始末端点，连成针线
  
  #判断针与直线是否相交
  if(trunc(NeedleY_1[i])!=trunc(NeedleY_2[i]))#判断针与直线不相交的逆事件
  {
    count_NeedleIntersect=count_NeedleIntersect+1
  }
  
  my_pi[i]=i/count_NeedleIntersect#计算累次针与直线相交之后π的近似值
}
my_pi=my_pi[!is.infinite(my_pi)]#去掉累计近似π的值的向量中的Inf值
piEst=times/count_NeedleIntersect#计算总共的针与直线相交之后π的近似值
piEst#输出π的近似值
plot(my_pi,pch=20,type ="b")#将所有可能的π的近似值画到图像上
abline(h=pi,col="red")#在散点图上画上y=π的直线

运行代码的结果一

以下为可视化结果，随着实验次数的增加，my_pi的值趋近于π。

R语言代码（优化版本）展示

以下给出线段与直线相交判定条件

根据对直线段半平面的定义，可给出如下引理：

引理        若任一直线段的两端点P、Q位于另一直线的同一半平面内，则线段不与直线相交，且必有。[1]

test<-function(LineSpacing,Needle_length,times,parallel_number)#LineSpacing为平行线间距，Needle_length为针的长度，times为投针次数，parallel_number为平行线的数目
{   if(LineSpacing>0&LineSpacing>Needle_length)
  {
      parallel_number=parallel_number-1#平行线的数量-1
      coord_Y<-parallel_number*LineSpacing#平行线的数目
      
      #以下为创建向量，即预分配内存，初始化向量
      my_pi=c()#初始化累计近似π的值的向量
      count<-0#初始化相交总数
      
      NeedleX_1=coord_Y*runif(times) #随机生成times次投针的起始点横坐标（乘上8为扩大随机0到1数为0到8）
      NeedleY_1=coord_Y*runif(times)#随机生成times次投针的起始点纵坐标（乘上8为扩大随机0到1数为0到8）
      NeedleLine_angle=pi*runif(times)#随机生成针与线的夹角
      NeedleX_2=NeedleX_1+cos(NeedleLine_angle)*Needle_length#计算投针的另一端点横坐标（三角函数）
      NeedleY_2=NeedleY_1+sin(NeedleLine_angle)*Needle_length#计算投针的另一端点纵坐标（三角函数）
      NeedleX<-c(NeedleX_1,NeedleX_2)#合并起始点横坐标和终点横坐标
      NeedleY<-c(NeedleY_1,NeedleY_2)#合并起始点纵坐标和终点纵坐标
      par(mfrow=c(1,2))#将绘图窗口划为1行两列，同时显示两图
      plot(NeedleX,NeedleY,pch=".")#绘制点图
      abline(h=seq(0,coord_Y,by=LineSpacing),col="red")#绘制平行线
      #连接始末端点，连成针线
      segments(NeedleX_1,NeedleY_1, NeedleX_2, NeedleY_2,col = "blue")
      #求针与直线相交的次数
      for (i in 1:times) 
      {
        for (j in seq(0,coord_Y,by=LineSpacing))
        {
          
          #每根针与每一条直线相交的判别（这里利用了线段与直线相交的判定条件）
          if((NeedleY_1[i]-j)*(NeedleY_2[i]-j)<=0)
          {
            count=count+1#每相交一次，总数就加一次
          }
        }
        my_pi[i]=(2*Needle_length*i)/(LineSpacing*count)#计算累次针与直线相交之后近似π的值
      }
      piEst=(2*Needle_length*times)/(LineSpacing*count)#计算总共的针与直线相交之后近似π的值
      piEst#输出近似π的值
      my_pi=my_pi[!is.infinite(my_pi)]#去掉累计近似π的值的向量中的Inf值
      plot(my_pi,pch=20,type ="b")#将所有可能的π的近似值画到图像上
      abline(h=pi,col="red")#在散点图上画上y=π的直线
      return(paste("近似π的值为",piEst))
  }
  else
    print("您输入的平行线间距小于零或者平行线间距小于针的长度,请重新输入test(平行线间距,针的长度,投针次数,平行线的数目)代码")
}
test(5,1,5000,10)#请输入数字，第一个为平行线间距，第二个为针的长度，第三个为投针次数，第四个为平行线的数目
test(4,2,4000,7)
test(3,1,10000,8)
test(5,3,10000,11)
test(8,3,10000,16)

 运行代码结果二

生成的数值和可视化图形（1）

> test(5,1,5000,10)
[1] "近似π的值为 3.19488817891374"

生成的数值和可视化图形（2） 

> test(4,2,4000,7)
[1] "近似π的值为 3.17208564631245"

> test(3,1,10000,8)
[1] "近似π的值为 3.31840053094409"

 

> test(5,3,10000,11)
[1] "近似π的值为 3.13234142521535"

 

> test(8,3,10000,16)
[1] "近似π的值为 3.05002033346889" 

 

参考文献

[1]夏德麟,熊鑫.线段相交性问题求解的新算法与原理[J].计算机辅助设计与图形学学报,1996,(02):95-100.