How to calculate the statistical uncertainty in a particle physics simulation?

66 Views Asked by At

I have a Monte Carlo code which simulates ions in a Tokamak. Most of the particles remain trapped forever. However, some of the particles escape. I can use my code to predict the fraction of particles that escape versus those that stay trapped. The code models collisions using a stochastic collision operator. Therefore, there is a statistical uncertainty with my results. Do you know any techniques to calculate this statistical uncertainty?

One idea I had is to use the Bootstrapping technique. I'll quickly explain my plan now:

Suppose we have $N$ trapped particles and $n<N$ that escape. We denote the particles that are trapped with $t_i$ (for $i\in\{1,2,...,N\}$) and the ones that escape with $e_i$ (for $i\in\{1,2,...,n\}$). Hence, the estimate for the fraction of markers that escape is $n/(n+N)$. To calculate the uncertainty, we make a list which looks like the following: $$\{t_1,t_2,...,t_N,e_1,e_2,...,e_n\}.$$ We resample this list (with replacement) $m$ times (e.g. $m=100$) to produce lists which may look something like this: $$\underbrace{\{t_{36}, e_{12}, t_{444}, t_{321},...\}}_{n+N\ \mathrm{elements}}.$$ Then calculate the fraction of markers which escape for each of the $m$ cases and take away $n/(n+N)$ to produce a list of values with length $m$. Finally, take, e.g. a 95% quantile of this to get a 95% confidence interval.

Can you see any problems with this method? Is there a more straightforward method I could be using? I know barely any statistics, so any advice is appreciated.