m = 10^5;  s = numeric(m);  rate=1/3;  lam=5
 for (i in 1:m) {
    n = rpois(1, lam)
    s[i] = sum(rexp(n+1, rate))  }
 mean(s);  var(s)
 ## 17.93764  # approximates E(S) = (5 + 1)(3)
 ## 97.76707  # approximates V(X)
 mean(s > 30)
 ## 0.11677   # approximates P(S > 30)