Sample size required for testing simulation results

Question

RE-WRITTEN QUESTION: I have a set of inequalities which forms a region in the 2D space. My main goal is to test whether the area of this region is larger or less than a threshold. Since the inequalities are difficult to solve by analytical method. I sample a set of points in the 2D space, and to justify how many points fall into that region formed by inequalities (it's easier to do so). Consequently, I can get an estimate of the region area. Now to satisfy accuracy (maybe a ratio of the region area, or an absolute value), I need to mathematically calculate out how many sample points in total I have to generate. So this is the main problem.

Now there is a whole space (whose area is 1), and a target region (whose area is $\theta$). I want to know whether the area of target region is larger than $\theta_0$ or less than $\theta_0$. Now I use hit-or-loss Monte Carlo method, which samples $N$ points from the whole space, and suppose we can get $M$ of them are in the target region. By which we obtain an estimate of $\theta$, i.e., $\hat{\theta} = \frac{M}{N}$.

Obviously, $\hat{\theta}$ is an unbiased estimate of $\theta$ with variance $\frac{\theta(1-\theta)}{N}$. Suppose given a confidence level $(1-\gamma)$, the estimate error satisfies $|\theta-\hat{\theta}|\leq z_{\gamma/2} \sqrt{\frac{\theta(1-\theta)}{N}} $. Now I want to determine a required sample size with an error bound $E$.

If the error bound is absolute error, i.e., $|\theta-\hat{\theta}|\leq z_{\gamma/2} \sqrt{\frac{\theta(1-\theta)}{N}} \leq E$, I can get \begin{align} N_1 \geq \frac{(z_{\gamma/2})^2 \theta (1-\theta)}{E^2}. \end{align}

And if the error bound is relative error, i.e., we need to make sure $\frac{|\theta - \hat{\theta}|}{\theta} \leq E$, then it equals to $|\theta-\hat{\theta}|\leq \theta E$, which means \begin{align} N_2 \geq \frac{(1-\theta)}{\theta}\frac{(z_{\gamma/2})^2}{E^2}. \end{align}

MY QUESTION IS: Since $\theta$ is unknown before, can I use $\theta_0$ (the threshold I want to judge) to substitute $\theta$ to get $N_1$ and $N_2$? Is this choice reasonable enough? What is the drawback of doing so? And what is the difference between using absolute error and relative error on the result of sample size?

Answer 1

It seems that you are using the standard Wald CI for an unknown binomial proportion $\theta$. This is an asymptotic CI, which is a poor choice for small $n.$ However, if you are doing a simulation with very large $n$ you should get reasonable results.

This style of CI assumes that $Z = \frac{\hat \theta - \theta}{\sqrt{\theta(1-\theta)/n}} \stackrel{aprx}{\sim} \mathsf{Norm}(0,1)$ (by the Central Limit Theorem) and then, because unknown $\theta$ is estimated by $\hat \theta = M/n$, that it is OK to use $\hat \theta$ for $\theta.$ Notice that you have two approximations: (a) assuming that $Z$ is normal and (b) assuming that $\hat \theta$ is near $\theta.$ Both assumptions are made more reasonable by large $n.$

If you want the worst-case scenario, you can estimate the margin of error using $\theta = 1/2$ because that gives the largest estimated margin of error (and so for you, the largest estimated sample size $n$ for a desired maximum margin of error).

If $\theta$ is very near $0,$ then the relative error can be very large. I have never tried getting a CI for relative error in such circumstances, but I suspect that for $\theta \approx 0$ there would be serious problems using $\hat \theta$ for $\theta.$ Notice the $\theta$ in the denominator of your last displayed equation.

---

Example: Simulating $\pi/4$ as the fraction of the unit circle within the square with vertices at $(0,1)$ and $(1,1)$ with R code. We simulate $B = 1,000,000$ points placed at random in the square, of which 785,161 within a unit of the origin are accepted, giving $\hat \theta = .7852$ compared with $\theta = \pi/4 = .7854.$ The absolute error is about .0002, which is well within the Wald 95% margin of error .0008.

B = 10^6;  u1 = runif(B);  u2 = runif(B)
acc = (u1^2+u2^2<1);  th.hat= mean(acc);  th.hat; pi/4;  abs(th.hat - pi/4)
## 0.785161                  # th.hat
## 0.7853982                 # th = pi/4
## 0.0002371634              # absolute error
1.96*sd(acc)/sqrt(B)      
## 0.0008049932              # approx margin of error from simulation
1.96*sqrt(th.hat*(1-th.hat)/B)
## 0.0008049928              # Wald margin of error

---

Addendum (responding to comments and revision): In theory, you are asking something impossible. But maybe there is a solution that will work in practice.

Take my simulation as an example. If you want to be 95% confident that $|\hat \theta - \theta| < \Delta$, then you want the 95% CI $\hat \theta \pm 1.96\sqrt{\hat\theta(1-\hat\theta)/N}$ to have $1.96\sqrt{\hat\theta(1-\hat\theta)/N} < \Delta.$ If $\Delta = 0.001$ and $\theta \approx 0.785$ then you can solve to get $N \approx 800,000$. And in my simulation with $N = B = 1,000,000$ I got a margin of error about 0.0008. As it turned out my simulated value $\hat \theta$ was within 0.0002 of $\theta_0 = \pi/4.$ However, knowing the exact value of $\pi$ is a famous mathematical quest. So if I had gotten $\hat \theta = 0.7853982,$ then I would have to know the value of $\pi/4$ to many places to know whether $\hat\theta > \pi/4$ or $\hat\theta < \pi/4.$

For planning purposes, you should pick a small $\Delta$ that seems good enough, and then find $N$ such $z_{\gamma/2}\sqrt{\theta_0(1-\theta_0)/N} < \Delta.$ So you can very likely come "within $\Delta$" of $\theta_0$ with your simulation. (For 'insurance' you might use $\theta_0 = 1/2$ to get the largest $N$ according to this method.) But from simulation and statistical analysis of simulation results, you cannot know for sure whether your simulated $\hat \theta$ is a little above or a little below the true $\theta_0.$ (Otherwise, some clever person could have used simulation to find $\pi$ to infinitely many digits.)