用极坐标方法模拟随机变量

Question

我有以下算法

步骤1.生成u1和u2~U（0,1）

步骤2.定义v1 = 2u1-1，v2 = 2u2-1和s = v1 ^ 2 + v2 ^ 2

步骤3.如果s> 1，请返回步骤1。

步骤4.如果s <= 1，则x = v1（-2logs / s）^（1/2），y = v2（-2logs / s）^（1/2）

以下是我在R中实现此算法的方法：

    PolarMethod1<-function(N)
{

  x<-numeric(N)
  y<-numeric(N)
  z<-numeric(N)

  i<-1

  while(i<=N)
  {u1<-runif(1)
  u2<-runif(1)
  v1<-(2*u1)-1
  v2<-(2*u2)-1
  s<-(v1^2)+(v2^2)

  if(s<=1)
  {
    x[i]<-((-2*log(s)/s)^(1/2))*v1
    y[i]<-((-2*log(s)/s)^(1/2))*v2
    z[i]<-(x[i]+y[i])/sqrt(2) #standarization
    i<-i+1
  }
  else
    i<-i-1
  }

  return(z)
}
z<-PolarMethod1(10000)
hist(z,freq=F,nclass=10,ylab="Density",col="purple",xlab=" z values")
curve(dnorm(x),from=-3,to=3,add=TRUE)

幸运的是，代码没有标记任何错误，并且与N=1000工作得很好但是当我改为N=10000时，而不是更好地处理曲线显示：

与N = 1000显示对比：

这是为什么？

我的代码有问题吗？当N增加时，它应该被更好地调整。

注意：我在代码中添加了z，以在输出中包含这两个变量。

Answer 1

为什么1000到100000次运行之间存在差异？

运行1000次模拟时，z值通常从-3.2到3.2。但是如果你将运行增加到100k你将获得更多的极值，z将从-4变为4。

直方图将z结果合并为10个区间。 z中较高的范围将导致较宽的箱，而较宽的箱通常会对概率密度进行调整。

1000次运行的箱宽约为0.5，但100k为1。

Answer 2

绘制直方图时要求10个箱，但这只是一个建议。你实际上得到了8，因为要覆盖从-4到4的范围，没有分成10个分区，最终得到漂亮的圆形数字，而8个分区有非常好的边界。

如果你想要更多的箱子，那么不要指定nclass。默认给了我20个箱子。或者指定breaks = "Scott"，它使用不同的规则来选择箱子。我看到使用此选项的大约80个箱子。