模拟复合随机变量 S答案

【问题标题】：Simulate the compound random variable S模拟复合随机变量 S
【发布时间】：2021-04-12 14:02:18
【问题描述】：

令 S=X_1+X_2+...+X_N 其中 N 是非负整数值随机变量，X_1,X_2,... 是 i.i.d 随机变量。（如果 N=0，我们设置 S=0）。

在 N ~ Poi(100) 和 X_i ~ Exp(0.5) 的情况下模拟 S。（绘制直方图并使用 numpy 或 scipy 内置函数）。并检查方程 E(S)=E(N)*E(X_1) 和 Var(S)=E(N)*Var(X_1)+E (X_1)^2 *Var(N)

我试图解决它，但我还不确定一切，并且还卡在了直方图部分。注意：我是 python 的新手，或者更一般地说，我是编程新手。

我的作品：

import scipy.stats as stats
import matplotlib as plt
N = stats.poisson(100)
X = stats.expon(0.5)
arr = X.rvs(N.rvs())
S = 0
for i in arr:
    S=S+i
print(arr)
print("S=",S)
expected_S = (N.mean())*(X.mean())
variance_S = (N.mean()*X.var()) + (X.mean()*X.mean()*N.var())
print("E(X)=",expected_S)
print("Var(S)=",variance_S)

【问题讨论】：

标签： python random simulation exponential poisson

【解决方案1】：

您现有的代码看起来很合理，但我会简化：

arr = X.rvs(N.rvs())
S = 0
for i in arr:
    S=S+i

到：

S = X.rvs(N.rvs()).sum()

要绘制直方图，您需要来自该分布的许多样本，现在可以通过以下方式轻松完成：

arr = []
for _ in range(10_000):
  arr.append(X.rvs(N.rvs()).sum())

或者，等效地，使用list comprehension：

arr = [X.rvs(N.rvs()).sum() for _ in range(10_000)]

要将这些绘制在直方图中，您需要 Matplotlib 中的 pyplot 模块，因此您的导入应该是：

from matplotlib.pyplot import plt

plt.hist(arr, 50)

上面的 50 表示在绘制直方图时使用该数量的“bin”。我们还可以将这些与您通过假设分布很好地近似于正态来计算的均值和方差进行比较：

approx = stats.norm(expected_S, np.sqrt(variance_S))

_, x, _ = plt.hist(arr, 50, density=True)
plt.plot(x, approx.pdf(x))

这是可行的，因为从matplotlib's hist method 返回的第二个值是垃圾箱的位置。我使用了density=True，因此我可以使用概率密度，但另一种选择可能是将密度乘以样本数以获得预期的计数，就像之前的直方图一样。

运行它给了我：

【讨论】：