在R中生成概率密度函数与计算期望及球面随机点距离统计问题

先给你理清楚数学逻辑，再结合你的sample3d函数做模拟验证：

先搞理论推导（省得模拟瞎跑） #

假设地球是半径为R的球体，球面上随机两点的球心角θ（两点和球心连线的夹角）的概率密度是：

f_θ(θ) = (1/2)sinθ，θ ∈ [0, π]

而你要的欧氏距离D（穿地球内部的直线距离）和θ的关系是：D = 2R sin(θ/2)。做变量替换后，可以推导出D的PDF：

f_D(D) = D/(2R²)，D ∈ [0, 2R]

对应的期望E(D)直接积分就能算出来：

E(D) = 4R/3

如果是单位球（R=1），那期望就是4/3≈1.333；地球半径取6371km的话，期望就是4*6371/3≈8494.67km。

用你的sample3d函数做模拟验证 #

假设你的sample3d函数是生成单位球面上的随机点（如果不是，最后调整R值就行），咱们一步步来：

步骤1：生成点对并计算D #

set.seed(456)
n <- 100000  # 样本量越大，结果越准
# 生成n个随机点
points <- sample3d(n)
# 再生成另一组独立的n个点
points2 <- sample3d(n)
# 计算每对点的欧氏距离：√[(x1-x2)²+(y1-y2)²+(z1-z2)²]
D <- sqrt(rowSums((points - points2)^2))

步骤2：估计D的PDF（和理论对比） #

# 用核密度估计PDF
kde_D <- density(D)
# 画图对比估计值和理论值（单位球理论PDF是D/2）
plot(kde_D, lwd=2, col="blue", main="PDF of Euclidean Distance D (Unit Sphere)",
     xlab="D", ylab="Probability Density")
# 添加理论PDF曲线
x_seq <- seq(0, 2, length.out=100)
y_theory <- x_seq/2
lines(x_seq, y_theory, lwd=2, col="red", lty=2)
# 加图例区分
legend("topright", legend=c("Kernel Density Estimate", "Theoretical PDF"),
       col=c("blue", "red"), lty=c(1,2), lwd=2)

步骤3：估计期望E(D) #

# 用样本均值估计期望
E_D_estimate <- mean(D)
cat("D的期望估计值：", E_D_estimate, "\n")
# 单位球的理论期望
cat("单位球的理论期望：", 4/3, "\n")

如果你的sample3d生成的是半径为R的球面上的点，只需要把理论公式里的R替换成你的球半径就行，模拟结果会自动贴合理论值。

内容的提问来源于stack exchange，提问作者user9519686