模拟连续随机变量的期望

Question

目前我想生成一些样本以获得它的期望和方差。

给定概率密度函数：f(x) = {2x, 0 <= x <= 1; 0否则}

我已经发现 E(X) = 2/3，Var(X) = 1/18，我的详细解法来自这里 https://math.stackexchange.com/questions/4430163/simulating-expectation-of-continuous-random-variable

但这是我在使用 python 进行模拟时所得到的：

import numpy as np
N = 100_000
X = np.random.uniform(size=N, low=0, high=1)
Y = [2*x for x in X]
np.mean(Y) # 1.00221 <- not equal to 2/3
np.var(Y) # 0.3323 <- not equal to 1/18

我在这里做错了什么？先谢谢你了。

Answer 1

为了逼近 x 的某些函数的积分，例如 g(x)，超过 S = [0, 1]，使用 Monte Carlo 模拟，你

在[0, 1]中生成N个随机数（即从均匀分布U[0, 1]中抽取）
计算 g(x_i) 在 i = 1 到 i = N 的算术平均值，其中 x_i 是第 i 个随机数：即(1 / N) 乘以 g(x_i) 的 i = 1 到 i = N 的总和。

第2步的结果是积分的近似值。

具有 pdf f(x) 和一组可能值 S 的连续随机变量 X 的期望值是 x * f(x) 在 S 上的积分。 X的方差是X的期望值-平方减去X的期望值的平方。

期望值：近似计算x * f(x)对S = [0, 1]的积分（即期望值 X)，设置 g(x) = x * f(x) 并应用上述方法。
方差：逼近(x * x) * f(x)对S = [0, 1]的积分（即期望值X-平方），设置 g(x) = (x * x) * f(x) 并应用上述方法。将此结果减去 X 的预期值估计值的平方，以获得 X.

方差

调整你的方法：

import numpy as np
N = 100_000
X = np.random.uniform(size = N, low = 0, high = 1)

Y = [x * (2 * x) for x in X]
E = [(x * x) * (2 * x) for x in X]

# mean
print((a := np.mean(Y)))
# variance 
print(np.mean(E) - a * a)

输出

0.6662016482614397
0.05554821798023696

与其制作 Y 和 E 列表，更好的方法是

Y = X * (2 * X)
E = (X * X) * (2 * X)

Y、E 在这种情况下是 numpy 数组。这种方法更有效。尝试制作 N = 100_000_000 并比较两种方法的执行时间。第二个应该快得多。

Answer 2

当您需要 X 本身的均值和方差时，您正在生成 Y = 2X 的均值和方差。你知道密度，但是 CDF is more useful for random variate generation than the PDF。对于你的问题，密度是：

所以 CDF 是：

鉴于 CDF 是范围 [0,1] 的易逆函数，您可以使用 inverse transform sampling 通过设置 F(X) = U 来生成 X 值，其中 U 是一个 Uniform( 0,1) 随机变量，并反转关系以求解 X。对于您的问题，这会产生 X = U^1/2.

换句话说，你可以用

生成X值

import numpy as np
N = 100_000
X = np.sqrt(np.random.uniform(size = N))

然后用数据做任何你想做的事，例如计算均值和方差、绘制直方图、在模拟模型中使用，或其他任何事情。

直方图将确认生成的数据具有所需的密度：

import matplotlib.pyplot as plt

plt.hist(X, bins = 100, density = True)
plt.show()

生产

然后可以直接根据数据计算均值和方差估计值：

print(np.mean(X), np.var(X))     # => 0.6661509538922444 0.05556962913014367

但是等等！还有更多...

误差范围

模拟生成随机数据，因此均值和方差的估计值在重复运行中会发生变化。统计学家使用 confidence intervals 来量化统计估计中不确定性的大小。当样本量大到足以调用中心极限定理时，均值的区间估计值计算为 (x-bar ± half-width)，其中 x-bar 是均值的估计值。对于 so-called 95% 的置信区间，half-width 是 1.96 * s / sqrt(n) 其中：

s是估计的标准差；
n是用于估计均值和标准差的样本数；和
1.96 是从正态分布和所需置信度派生的比例常数。

half-width 是误差范围的定量度量，a.k.a。精度，估计值。请注意，随着 n 变大，估计的误差范围更小并且变得更精确，但是由于平方根，样本量的增加会减少 returns。如果使用独立抽样，将精度提高 2 倍将需要 4 倍的样本量。

在Python中：

var = np.var(X)
print(np.mean(X), var, 1.96 * np.sqrt(var / N))

产生的结果如

0.6666763186360812 0.05511848269208021 0.0014551397290634852

其中第三列是置信区间 half-width。

提高精度

如果我们使用基于期望和方差的基本属性的巧妙技巧，逆变换采样可以为给定的样本量产生更高的精度。在介绍 prob/stats 课程中，您可能被告知 Var(X + Y) = Var(X) + Var(Y)。真正的关系实际上是 Var(X + Y) = Var(X) + Var(Y) + 2Cov(X,Y)，其中 Cov(X,Y) 是 X 和 Y 之间的 covariance。如果它们是独立的，则协方差为 0，一般关系成为我们 learn/teach 在介绍课程中的关系，但如果它们不是独立的，则必须使用更一般的方程。方差总是一个正数，但协方差可以是正的也可以是负的。因此，很容易看出，如果 X 和 Y 具有负协方差，则它们之和的方差将小于它们独立时的方差。负协方差意味着当 X 高于其均值时 Y 倾向于低于其均值，而 vice-versa.

那么这有什么帮助呢？这很有帮助，因为我们可以使用逆变换以及一种称为 antithetic variates 的技术来创建成对的随机变量，这些变量分布相同但具有负协方差。如果 U 是具有 Uniform(0,1) 分布的随机变量，则 U' = 1 - U 也具有 Uniform(0,1 ）分配。（事实上，翻转任何对称分布都会产生相同的分布。）因此，X = F^-1(U) 和 X' = F^-1(U') 同分布，因为它们由相同的 CDF 定义，但具有负协方差，因为它们落在他们共享中位数的相反两侧，因此强烈倾向于落在他们的平均值的相反两侧。如果我们对每对进行平均得到 A = (F^-1(u_i) + F^{- 1}(1-u_i)) / 2)期望值E[A] = E [(X + X')/2] = 2E[X]/2 = E[X] 而方差 Var(A) = [(Var(X) + Var(X') + 2Cov( X,X')]/4 = 2[Var(X) + Cov(X,X')]/ 4 = [Var(X) + Cov(X,X')]/2。换句话说，我们得到一个随机变量A 其平均值是 X 均值的无偏估计，但方差较小。

为了公平地比较对立结果 head-to-head 与独立抽样，我们采用原始样本大小并将其分配给 U[=122= 的逆变换生成的一半数据]'s，另一半是使用 1-U's 通过对偶配对生成的。然后我们像以前一样擦除配对值并生成统计信息。在 Python:

U = np.random.uniform(size = N // 2) antithetic_avg = (np.sqrt(U) + np.sqrt(1.0 - U)) / 2 anti_var = np.var(antithetic_avg) print(np.mean(antithetic_avg), anti_var, 1.96*np.sqrt(anti_var / (N / 2)))

产生的结果如

0.6667222935263972 0.0018911848781598295 0.0003811869837216061

请注意，使用独立采样产生的 half-width 几乎是使用对立变量产生的 half-width 的 4 倍。换句话说，我们需要多一个数量级的独立采样数据才能达到相同的精度。

模拟连续随机变量的期望

Simulating expectation of continuous random variable

python

numpy

statistics

simulation

distribution

但是等等！还有更多...

误差范围

提高精度