KAHIBARO
22.5. Statistical Analysis
Table of Contents
Mean
In GATE data analysis, you often summarize large sets of numbers such as deposited energy values, hit positions, or timing information. The mean is the simplest statistical summary and represents the central or average value of a quantity.
Suppose you have a set of $N$ values, for example the energy of $N$ detected events, written as $x_1, x_2, \dots, x_N$. The arithmetic mean is defined as
$$
\bar{x} = \frac{1}{N} \sum_{i=1}^{N} x_i.
$$
In Python, once you have loaded data from GATE output into an array (for instance using NumPy), the mean is usually computed with a single function call, but it is important to understand what it represents. The mean of an energy spectrum gives you the average detected energy, the mean of a dose distribution in a region gives the average dose to that region, and the mean of a timing distribution can indicate a central detection time or a central time difference for coincidence events.
For Monte Carlo simulations, the mean of a quantity such as energy deposition per event is often used as an estimator of the true physical expectation value. When you repeat a simulation with independent random seeds, the different sample means you obtain will fluctuate around the true mean. As the number of events $N$ increases, these fluctuations become smaller and the sample mean becomes a more precise estimate.
It is sometimes useful to distinguish between an unweighted mean and a weighted mean. In many simple GATE analyses, each event has the same statistical weight and you use the unweighted mean defined above. In more advanced cases, such as phase space reuse or specialized variance reduction techniques, each event can carry a weight $w_i$, and then the weighted mean is
$$
\bar{x}_w = \frac{\sum_{i=1}^{N} w_i x_i}{\sum_{i=1}^{N} w_i}.
$$
Weighted means appear when one event in the simulated data represents a different importance compared to others. In standard introductory GATE simulations without variance reduction, you usually work with unweighted events and can stay with the simple definition.
Key idea: The mean, $\bar{x} = \frac{1}{N}\sum_{i=1}^{N} x_i$, is your basic estimator of the central value of a simulated quantity such as energy, dose, or time. As $N$ grows, the mean becomes a more reliable estimate of the true physical expectation value.
Standard deviation
While the mean describes the central value, the standard deviation describes how much the values fluctuate around that mean. In Monte Carlo simulations, this spread reflects both the physical variability of the process and the statistical uncertainty due to a finite number of events.
Given $N$ values $x_1, x_2, \dots, x_N$ with mean $\bar{x}$, the sample standard deviation $s$ is defined as
$$
s = \sqrt{\frac{1}{N - 1} \sum_{i=1}^{N} (x_i - \bar{x})^2}.
$$
The standard deviation has the same units as the original quantity. If $x_i$ are energies in keV, then $s$ is also in keV. A narrow energy distribution has a small standard deviation, while a broad distribution has a large one. The same applies to dose values, timing values, or any other quantity you analyze.
The standard deviation tells you about the spread of individual measurements or events, but when you run a GATE simulation you are usually interested in how precisely you know the mean. For that purpose, you often compute the standard error of the mean, which is
$$
\sigma_{\bar{x}} = \frac{s}{\sqrt{N}}.
$$
This quantity decreases when you use more events, since the uncertainty on the mean scales like $1 / \sqrt{N}$. In dose calculations, for example, you might compute for each voxel the mean dose and an associated uncertainty, often estimated using this formula on the energy deposition per event in that voxel.
A common way to summarize the uncertainty in Monte Carlo results is to quote the mean value together with one standard deviation of the mean. For example, you might report a voxel dose as $1.25 \pm 0.03$ Gy. Here 1.25 Gy is the estimated mean dose and 0.03 Gy is the estimated statistical uncertainty from a finite number of simulated histories.
In practice, you will usually compute standard deviations and associated uncertainties using numerical libraries such as NumPy, but it is important to know that these quantities are derived directly from the fluctuations of event level data that come from GATE.
Important formulas:
Standard deviation:
$$
s = \sqrt{\frac{1}{N - 1} \sum_{i=1}^{N} (x_i - \bar{x})^2}
$$
Standard error of the mean:
$$
\sigma_{\bar{x}} = \frac{s}{\sqrt{N}}
$$
These are central to quantifying Monte Carlo statistical uncertainty in GATE analysis.
Event counts
In GATE simulations, many basic statistical quantities are built from simple event counts. An event count is simply how many times a certain type of outcome occurs. Typical examples include the number of detected singles above an energy threshold, the number of coincidences in a PET simulation, the number of photons transmitted through a shield, or the number of particle tracks entering a region.
Let $N$ be the total number of primary events in the simulation. If you are interested in how many times a particular condition $C$ is satisfied, such as energy in a certain window, you count how many events satisfy it. Suppose $n$ events satisfy $C$. You can then form a fraction or probability estimate
$$
p = \frac{n}{N}.
$$
Counts and such ratios appear everywhere in GATE based analyses. For example, in PET you may compute the number of true coincidences, the number of scattered coincidences, and the number of random coincidences. By dividing each of these by the total number of coincidences, you obtain relative fractions such as scatter fraction. In shielding studies, you count how many photons are transmitted through a material compared with how many are incident on it to estimate transmission and attenuation.
Because these are Monte Carlo estimates, they have statistical uncertainties. For a simple count of $n$ events out of $N$ with probability $p$, the standard deviation of the count due to Poisson statistics is approximately
$$
\sigma_n \approx \sqrt{n},
$$
and the standard deviation of the fraction $p = n / N$ is approximately
$$
\sigma_p \approx \sqrt{\frac{p (1 - p)}{N}}.
$$
When the counts are large, these approximations are accurate and very useful. They tell you that if you want to reduce the relative uncertainty on a count by a factor of 2, you must increase the number of events by a factor of 4, because the relative uncertainty behaves roughly like $1 / \sqrt{n}$.
Many performance metrics used in medical physics simulation, such as sensitivity, detection efficiency, or transmission factors, are built from event counts and their ratios. The precision of these metrics is controlled by the number of events you simulate and the way you select which events to include in each count.
In practice, you obtain event counts by applying logical conditions in your analysis code to select subsets of data from the GATE output. You then summarize these selections by reporting both the counts and their associated statistical uncertainties. This provides a quantitative understanding of how reliable your simulated estimates are.
Key counting rules:
For a count $n$ of events:
Standard deviation of the count (Poisson):
$$
\sigma_n \approx \sqrt{n}
$$
For a fraction $p = n / N$:
Standard deviation of the fraction:
$$
\sigma_p \approx \sqrt{\frac{p (1 - p)}{N}}
$$
Event counts and these uncertainties are the foundation of most performance and efficiency metrics in GATE simulations.
Views: 11