# sample size Monte Carlo experiment

I'm studying Monte Carlo analysis but I find very counter-intituive the computation of the minimum sample size in order to reach a certain level of precision. As stated in Monte Carlo methods in Financial Engineering, such sample size comes from the CLT and it can be computed according to

$$n = \frac{z^2_{\delta/2} \sigma^2}{\epsilon^2}$$

My experiment is the classic Bernoulli variable, with $p$ probability of success. Suppose my level of precision is $0.01$, then my $z_{0.01/2} = 2.57$.

I consider three cases.

$p=0.1\%$

$p=1\%$

$p=10\%$

With $\sigma^2 = p(1-p)$, I get

$n_{10\%} = 2.57^2 \cdot 0.1 \cdot 0.9 / 0.01^2 = 5944$

$n_{1\%} = 2.57^2 \cdot 0.01 \cdot 0.99 / 0.01^2 = 654$

$n_{0.1\%} = 2.57^2 \cdot 0.001 \cdot 0.999 / 0.01^2 = 66$

Intuitevely I would think the rarer the event, the greater the number of the simulations needed in order to approach the theoretical mean. What am I getting wrong?

tl;dr: you're wrong, and your error is that you're mistaking absolute error and relative error.

Long response: the formula you reports gives the sample size needed to get a certain absolute estimation error. If the probability becomes smaller, it's obvious that you need less runs to get an estimate with the same absolute error, because intuitively your estimate will usually be close to zero, and thus to the quantity you're estimating. If I want to estimate a probability $p=0.01$ with an absolute error $\epsilon = 0.01$, i.e., a relative error of $100\%$, this means that 0 or 0.02 are acceptable estimates. Intuitively, this is much simpler than estimating a larger probability $p=0.5$ with the same absolute error $\epsilon = 0.01$, which in this case corresponds to a relative error of $2\%$.

Practical example: suppose $p=0.01$ is the probability of obtaining "heads" in a loaded coin, and I throw my coin just once. I will likely get "tails", since $p$ is so small. Now, my Monte Carlo estimate of $p$ is

$$\hat{p}_N=\frac{\sum_{i=1}^N X_i}{N}=\frac{\sum_{i=1}^1 0}{1} = 0$$

Is this acceptable, according to my absolute error criterion of $\epsilon = 0.01$? Obviously yes:

$$|\hat{p}_1-p|=|0-0.01|=0.01=\epsilon$$

Thus, as you can see, with just 1 throw of my coin, I got an acceptable (according to your criterion) estimate of $p$. Since the probability of getting "heads" in one throw is $0.01$, you can easily convince yourself that you will get an acceptable estimate on average $99$ times out of $100$, i.e., the confidence level (what you called precision) is $0.01$. On the other hand, note that the relative error is huge: in the best case, my one-throw estimate gives me a relative error

$$\frac{|\hat{p}_1-p|}{p}=\frac{0.01}{0.01}=1> 0.01=\epsilon$$

Consider instead the case of $p=0.5$: if you throw your coin only once, you can get an estimate of $0$ or $1$ with the same probability. This means that, on a single throw, your absolute error will be much larger than $0.01$ (it will always be $0.5$): you need more throws to get the desired absolute error.

You would have been right if you had considered the relative error: in that case, it's true that to get the same relative error when $p\to0$, you need much more throws than when $p$ is large. As a matter of fact, you can easily prove, following the same asymptotic approximation used to derive the formula you report, that the sample size needed to estimate the success probability of a Bernoulli RV with a relative estimation error $\epsilon$ and a confidence $\delta$ is

$$n=\frac{z^2_{\delta/2} \frac{(1-p)}{p}}{\epsilon^2}$$

I leave it as an exercise for you to prove that the paradox you observed before, has now disappeared.

However, note that in practice this formula is a very bad estimate of the sample size, because the CLT approximation for the sample mean of a sequence of Bernoulli RVs breaks down when $p\to 0$. Or, more precisely, the sample size needed in order for the normal approximation to be an accurate approximation of the distribution of $\hat{p}_N$, grows indefinitely when $p\to0$.

UPDATE: so, if the confidence interval for $\hat{p}_N$, based on the CLT, is so bad when $p\to 0$, what are better solutions? I include here some suggestions which I and the OP originally discussed in the comments, since comments should be temporary and only have the purpose to improve Q&A.

1. In the regime of large $n$ and small $p$, when the normal approximation to the binomial distribution becomes increasingly worse, the Poisson approximation becomes increasingly more accurate, so that could be a tool. See here

2. However, using another analytical approximation and then having to derive the corresponding CI by hand is surely educational, but hardly a good way of spending limited time if your goal is to quickly have an accurate estimate. Also, large $n$ and small $p$ or large $n$ and $p\approx 0.5$ are just two "boundary" regimes: what about other cases? Which CI should you use there? Of course, a lot of research has been devoted to the definition of CIs for binomial proportions which are more accurate and robust than the Wald (CLT-based) one. See Brown et al., Interval Estimation for a Binomial Proportion, Statistical Science, Volume 16, Issue 2 (2001), 101-133 for the related theory and

https://stackoverflow.com/questions/21719578/confidence-interval-for-binomial-data-in-r

Confidence interval for Bernoulli sampling

Binomial confidence interval estimation - why is it not symmetric?

to compute such intervals in R.

3. The original question, however, wasn't how to compute accurate & robust CIs for binomial proportions, or it would have been a duplicate. It was about computing the sample size needed for a certain relative or absolute accuracy $\epsilon$ in the estimation of $\hat{p}_N$. Thus, to get such an estimate directly without even having to invert the confidence interval, you can use this R package.

• I got where I was mistaken, thanks...about your last remark, could I use the Poisson approximation in order to get better results? – Vittorio Apicella Oct 3 '17 at 6:24
• Really thanks for the links. I've also found an R module which directly computes the correspondent sample size for many of the alternative estimators, it's called "binomSamSize". cran.r-project.org/web/packages/binomSamSize/index.html (hope it'll be useful also for someone else) – Vittorio Apicella Oct 3 '17 at 21:55