Search code examples
pythonrandomsimulationexponentialpoisson

Simulate the compound random variable S

Let S=X_1+X_2+...+X_N where N is a nonnegative integer-valued random variable and X_1,X_2,... are i.i.d random variables.(If N=0, we set S=0).

Simulate S in the case where N ~ Poi(100) and X_i ~ Exp(0.5). (draw histograms and use the numpy or scipy built-in functions).And check the equations E(S)=E(N)*E(X_1) and Var(S)=E(N)*Var(X_1)+E(X_1)^2 *Var(N)

I was trying to solve it, but I'm not sure yet of everything and also got stuck on the histogram part. Note: I'm new to python or more generally , new to programming.

My work:

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

  • Your existing code mostly looks sensible, but I'd simplify:

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

    down to:

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

    To draw a histogram, you need many samples from this distribution, which is now easily accomplished via:

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

    or, equivalently, using a list comprehension:

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

    to plot these in a histogram, you need the pyplot module from Matplotlib, so your import should be:

    from matplotlib.pyplot import plt

plt.hist(arr, 50)

    The 50 above says to use that number of "bins" when drawing the histogram. We can also compare these to the mean and variance you calculated by assuming the distribution is well approximated by a normal:

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

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

    This works because the second value returned from matplotlib's hist method are the locations of the bins. I used density=True so I could work with probability densities, but another option could be to just multiply the densities by the number of samples to get expected counts like the previous histogram.

    Running this gives me:

    money shot