如何在r中生成概率密度函数和期望？

Question

任务：

苍蝇埃里克有一个朋友厄尼。假设两只苍蝇坐在独立的位置，均匀分布在地球表面。让 D 表示 Eric 和 Ernie 之间的欧几里得距离（即，在穿过地球内部的直线上）。

对 D 的概率密度函数做一个猜想，并给出一个 估计其期望值 E(D)。

到目前为止，我已经制作了一个在地球表面上生成两个点的函数，但我不确定下一步该做什么：

sample3d <- function(2)
  {
  df <- data.frame()
  while(n > 0){
    x <- runif(1,-1,1)
    y <- runif(1,-1,1)
    z <- runif(1,-1,1)
    r <- x^2 + y^2 + z^2
    if (r < 1){
      u <- sqrt(x^2+y^2+z^2)
      vector = data.frame(x = x/u,y = y/u, z = z/u)
      df <- rbind(vector,df)
      n = n- 1
    }
  }
  df
}
E <- sample3d(2)

Answer 1

这是一个有趣的问题。我将概述一种计算方法；我会把数学留给你。

1. 首先，我们修复了一个随机种子以实现可重复性。

   set.seed(2018);
   

2. 我们从单位球面上采样10^4点。

   sample3d <- function(n = 100) {
     df <- data.frame();
     while(n > 0) {
       x <- runif(1,-1,1)
       y <- runif(1,-1,1)
       z <- runif(1,-1,1)
       r <- x^2 + y^2 + z^2
       if (r < 1) {
         u <- sqrt(x^2 + y^2 + z^2)
         vector = data.frame(x = x/u,y = y/u, z = z/u)
         df <- rbind(vector,df)
         n = n- 1
       }
     }
     df
   }
   df <- sample3d(10^4);
   

   请注意，sample3d 效率不高，但这是一个不同的问题。

3. 我们现在从df 中随机抽取2 个点，计算这两个点之间的欧几里得距离（使用dist），然后重复这个过程N = 10^4 次。

   # Sample 2 points randomly from df, repeat N times
   N <- 10^4;
   dist <- replicate(N, dist(df[sample(1:nrow(df), 2), ]));
   

   正如@JosephWood 所指出的，数字N = 10^4 有点随意。我们使用bootstrap 来推导经验分布。对于N -> infinity，可以证明经验引导分布与（未知）总体分布相同（引导定理）。经验分布和总体分布之间的误差项为1/sqrt(N)，因此N = 10^4 应该导致1% 左右的误差。

4. 我们可以将得到的概率分布绘制为直方图：

   # Let's plot the distribution
   ggplot(data.frame(x = dist), aes(x)) + geom_histogram(bins = 50);
   

4. 最后，我们可以得到平均值和中位数的经验估计。

   # Mean
   mean(dist);
   #[1] 1.333021
   
   # Median
   median(dist);
   #[1] 1.41602
   

   这些值接近理论值：

   mean.th = 4/3
   median.th = sqrt(2)