How to calculate this in R and in infinite series for exact value?

Given a random variable X with mean , the expectation E(f(X)) is sometimes used to characterize the asymmetry of the distribution. In the case of the negative binomial with parameters r, p this becomes which is an infinite series, with expanded form f(0)p" + f(1)rap" + f(2) ," 2 - 9 ? p" + f (3 ) " (7 1) (7 2 ) 16 2q'p t ... This would be challenging to calculate directly, However, we can use a Monte Carlo method to estimate as follows: Assuming the law of large numbers holds for f(X), then if we generate a large enough sample 1, $2, T3, . . . I, from the negative binomial distribution, then the sample mean f(21) + f(12) + ...+ f (In) n (a sequence average) will be close to the expected value E(f(X)). Let's check Monte Carlo simulation in R: As an example let's assume we want to calculate E(X) where X has negative binomial distribution with parameters r = 4 and p = .61. Infinite series calculations show that the exact value is rq(1 + rq)/p2 = 3.9936/0.3721 - 10.73259876. But doing it this way is a challenge. Lets check Monte Carlo simulation in R. | utilize this code: which is particular for estimating E(X) the case that r = 4 and p = .61. > w x X [1] 4 1 1 36 4 9 16 16 36 16 > mean(x) [1] 13.9 > x mean(x) [1] 10.74023 At first I did only 10 repetitions so that you can see all individual samples. Then "mean" will generate the average. Note that the mean of 10 samples is not accurate enough. So I did 1000000 repetitions so that the LLN has "time" to take effect. Note that our estimate 10.74023 is extremely close to the exact value. Your problem: Given an negative binomial random variable X with parameters r = 16 and p = 0.21, use the Monte Carlo method to compute E((X -() ), where u is, of course, the mean of X given by u = rq/p. ( It might be easiest if you just modify the code above, using / either as an exact fraction or as a decimal rounded to 6 or more decimal places.) You will be marked correct is your answer is within 1% of the true value. (You might get it first try, but you might be unlucky and not. If not, you may need to run the code a few times to "resample".)