模拟复合随机变量 S

Question

令 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)

Answer 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 表示在绘制直方图时使用该数量的“bins”。我们还可以将这些与您计算的均值和方差进行比较，假设分布非常接近正态分布：

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 返回的第二个值是 bin 的位置。我使用 density=True 所以我可以使用概率密度，但另一种选择可能是将密度乘以样本数以获得预期的计数，如之前的直方图。

运行这给了我：

模拟复合随机变量 S

Simulate the compound random variable S

python

random

simulation

exponential

poisson