Counting statistics (Radioactivity) Widget

Counting statistics (Radioactivity) Widget

The Counting Statistics (Radioactivity) widget is an absolute goldmine of learning resource for anyone involved in the measurement of radiation / radioactivity at anything other than the most basic level (or just interested in the subject). Chris has written an extensive paper to explain counting statistics, and how the widget can be used to illustrate the concepts in an interactive way (without a real radioactive source).

Counting statistics

Dr Chris Robbins, Grallator

Monitoring radioactivity - background

The basis of all radioactivity monitoring essentially boils down to counting decay events. Radioactive decay is a random process and is described by Poisson statistics. The Poisson distribution is a discrete probability distribution that gives the probability of a given number of events occurring in a given interval of time. For a process that has an expected number of events, \(\lambda\), in a given time interval, the probability of exactly \(N\) events happening is given by \[ P(N)=\frac{\lambda^N e^{-\lambda}}{N!} \] Note, \(N\) has to be an integer; this is a discrete distribution - a decay either has or has not occurred (the detector acts as an observer - there is no Schrödinger's cat scenario!). However, the expected number of events \(\lambda\) does not have to be integer. This is obvious when you think about the definition of it being "...in a given interval of time"; an expected [integer] number of 975 in a 1 hour period translates to an expected number of 16.25 [not an integer] in a 1 minute period.

The figure below shows the general shape of the Poisson distribution for three values of \(\lambda\). The dotted lines are for clarity - the markers are the actual "output".

The key features are:

• The distribution is skewed, although the extent of the skew reduces for larger values of \(\lambda\).
• The mean value of the distribution is \(\lambda\).
• The standard deviation of the distribution is \(\sigma = \sqrt{\lambda}\).
• The mode value(s) are \( \lceil \lambda \rceil - 1, \lfloor \lambda \rfloor \), where \( \lceil \cdot \rceil\) denotes next integer above the value and \( \lfloor \cdot \rfloor\) denotes next integer below the value.

Counting events in a single count

In monitoring using a single count over a period of time it is assumed that the resulting count value, \(X\), is the equal to the actual mean value, \(\lambda\). This is justified as \(\lambda\) is (one of) the mode value(s), and so is a more likely outcome. As we are assuming Poisson statistics, the standard deviation of a count \(X\) is \(\sigma = \sqrt{X}\), and this is used to quantify the error in the count to allow for experimental statistical variation of \(X\). The actual count would therefore be reported as \[ X \pm \sqrt{X} \]

Aside - Poisson and normal

For the case of a normal distribution, the above would state that 68% of all measured counts would be in the range \(X \pm \sigma \), and wider confidence intervals could be constructed, e.g. the 95% confidence interval would be \( X \pm 1.96\sigma \). The Poisson distribution is not the normal distribution and such statements are not strictly true for the Poisson distributions. However, for large values of \(X\), a normal distribution approximates a Poisson distribution and the 68% and 95% confidence intervals can be approximated using \(\pm \sqrt{X} \) and \(\pm 1.96\sqrt{X} \). A "large" value of \(X\) is somewhat subjective; some use \(X \ge 20\), others \(X \ge 30\), and others large values still. (Note also, there are methods for calculating confidence intervals for small counts but these are beyond the scope of this text.)

The key take away point is that if you are going to use normal approximation confidence intervals, the monitoring period should be long enough that a large number of counts are made.

From count to count rate

A count value of \(X\) obtained over a time period \(t\) gives a count rate of \(R =\frac{X}{t} \) with a corresponding standard deviation of \( \sigma_R = \frac{\sqrt{X}}{t} \). The rate result is therefore reported as \[ R \pm \sigma_R \]

Using \(R =\frac{X}{t} \Leftrightarrow X = Rt \) in the expression for \( \sigma_R\) gives \[ \sigma_R = \frac{\sqrt{X}}{t} = \frac{\sqrt{Rt}}{t} = \sqrt{\frac{Rt}{t^2}} = \sqrt{\frac{R}{t}} \]

Example 1

A count of \(X=500\) is made over a time period of \(t=10\) minutes.

The rate is \(R =\frac{X}{t} =\frac{500}{10} = 50.000 \) counts per minute (cpm).

The error (one standard deviation) is \(\sigma_R = \sqrt{\frac{R}{t}} = \sqrt{\frac{50}{10}} = 2.236\).

The reported count rate is \(50.000 \pm 2.236\) cpm.

The sampling simulator can be used to recreate this sampling example, see image below.

In this set up the background rate and sampling time are set at zero (more on background later) and the source actual activity is set to 50 cpm. As the sliders are moved a new sample will automatically be taken. Repeated samples can be made with the settings unchanged by pressing the "Simulate sampling" button.

The particular sample in the above has recorded a value of \(51.400 \pm 2.267\). As might be expected the value isn't 50 but 50 is within the stated error range. Repeated resampling will give different count values and reasonably similar errors.

As the monitoring of radiation is measuring a random process there are times when a single measurement may yield unexpected results. For example, the image below shows a sampled count of \(46.400 \pm 2.154\).

As we are in control of the situation we know the actual rate is 50 but this is outside of the quoted error range on the value of 46.4, so what's happened? Recall the result is presented as \(X \pm \sigma \), i.e. the error is noted as one standard deviation. Approximating a Poisson distribution by a normal distribution (careful!) this would mean that 68% of all measurements would fall within the stated range, or 32% won't. This is a case that didn't! If you increase the monitoring time you will reduce the error - this won't prevent actual values being outside the sample error range when quoted at one standard deviation, but it will mean that if it is, the reported value will be closer to the actual value. (Another way of reducing the error is to take more samples, but this isn't discussed here.)

The image below shows monitoring for 120 minutes. Again, the actual value is outside the sampled value's [one standard deviation] error range but it's closer to the actual value than the previous case.

The cautionary note here is to be aware of what the \(\pm \space error\) means when a value is quoted. It doesn't necessarily ensure the actual value is within the range, it only states the actual value is within the stated range to within a given confidence level. If you want more confidence that the actual value is within the range you need to use multiples of the standard deviation, which necessarily leads to a larger reported error.

Obtaining a specified error level

The above has reported a rate value and assigned an error based on one standard deviation. The relative size of this error is given by \[ \% \sigma_R = 100 \times \frac{ \sigma_R}{R} \] For the values in Example 1: \[ \% \sigma_R = 100 \times \frac{ \sigma_R}{R} = 100 \times \frac{ 2.236}{50} = 4.472\% \]

The obvious question to ask is "I want the achieve an error of at least \(x\% \), how do I achieve this?". As you cannot control the rate of decay the only variable you can change to reduce the error is to increase the counting time. The time required for a specified value of \(\% \sigma_R\) is calculated as follows: \[ \% \sigma_R = 100 \times \frac{ \sigma_R}{R} \Rightarrow \sigma_R = \frac{\% \sigma_R R}{100} \] Recall, \(\sigma_R = \sqrt{\frac{R}{t}}\), so that \[ \begin{align} \sigma_R = \sqrt{\frac{R}{t}} &=\frac{\% \sigma_R R}{100} \\ \\ \frac{R}{t} &= {\left( \frac{\% \sigma_R R}{100} \right)}^2 \\ \\ t&=\frac{{100}^2 R}{\%\sigma_R^2 R^2} \\ \\ t&= \frac{1}{R} {\left( \frac{100}{\% \sigma_R} \right)}^2 \end{align} \]

For the values in Example 1, the required monitoring time to achieve a 1% error is \[ t = \frac{1}{R} {\left( \frac{100}{\% \sigma_R} \right)}^2 = \frac{1}{50} {\left( \frac{100}{1} \right)}^2 = 200 \space minutes \] Note how a reduction in the error by a factor 4.472 requires an increase in measurement time required by a factor of 20 ( \(=4.472^2\) ).

The above can be demonstrated using the simulator. Selecting a time value of 180 minutes gives a relative error of just over 1%, while using 240 minutes gives an error of just under 1%. An example with a 240 minute time is shown below.

Accounting for background radiation

Real world monitoring is complicated by the fact that the environment ususally has an intrinsic backgorund source of radiation that adds to the count of any source of interest. To determine the activity of just the source the background activity must be removed. Denfining the following:

• \(R_B\) as the measured background count rate.
• \(R_M\) as measured count rate from both the source and background radiation.
• \(R_S\) as the count rate associated with just the source (what we want to find).

The value of \(R_S\) is given by subtracting the background rate from the measured rate: \[ R_S = R_M - R_B \] The values of \(R_B\) and \(R_M\) will have associated uncertainies, \(\sigma_{R_B} \) and \(\sigma_{R_M} \). As the two measurements of \(R_B\) and \(R_M\) are independent of each other, the combined uncertainty is given by summing the individual uncertainties in quadratutre, so that \[ \sigma_{R_S} = \sqrt{\sigma_{R_B}^2 + \sigma_{R_M}^2} \]

Example 2

In the absence of a source, a count of \(X_B=300\) is made over a time period of \(t_B=10\) minutes, giving a rate \(R_B =\frac{X_B}{t_B} =\frac{300}{10} = 30.000 \) cpm.

The error (one standard deviation) on \(R_B\) is \(\sigma_{R_B} = \sqrt{\frac{R_B}{t_B}} = \sqrt{\frac{30}{10}} = 1.732\).

In the presence of a source, a count of \(X_B=500\) is made over a time period of \(t_M=10\) minutes, giving a rate \(R_M =\frac{X_M}{t_M} =\frac{500}{10} = 50.000 \) cpm.

The error (one standard deviation) on \(R_M\) is \(\sigma_{R_M} = \sqrt{\frac{R_M}{t_M}} = \sqrt{\frac{50}{10}} = 2.236\).

The net source only count rate is \(R_S = R_M - R_B = 50.000 - 30.000 = 20.000 \) cpm.

The source only error is \(\sigma_{R_S} = \sqrt{\sigma_{R_B}^2 + \sigma_{R_M}^2}=\sqrt{{1.732}^2 + {2.236}^2} = 2.828 \)

The reported source only rate is reported as \( R_S = 20.000 \pm 2.828\) cpm.

The above situation is reconstructed in the simulator by selecting a background rate of 30 and a source rate of 20. There are now two sampled results shown, one for a background only measurement and a second for a background plus source measurement. The constructed source only value, with associated error, is shown under the two sampled values.

Vary the values and see how the results change.

Obtaining a specified error level

As with Example 1, the question of how long should the measuring period be to obtain a specified error level can be asked. The answer isn't as straightforward as before as there are now two times to consider, \(t_B\) and \(t_M\). To answer the question of the total time, \(T\), required to obtain a given accuracy, where \(T=t_B+t_M\), we first need to consider the problem of how to optimally apportion the values \(t_B\) and \(t_M\) such that the error in the net source only value is as small as possible.

The optimal ratio of measurement times

It is required to choose \(t_B\) and \(t_M\) such that \(T=t_B+t_M\) and which minimises the source only error, \(\sigma_{R_S} = \sqrt{\sigma_{R_B}^2 + \sigma_{R_M}^2}\). As \(\sigma_{R_S}\) is, by definition, a positive value, minimising \(\varepsilon = \sigma_{R_S}^2\) will be equivalent to minimising \(\sigma_{R_S}\). The mimimum is found by taking the derivative of \(\varepsilon\) and setting to zero (and using \(\sigma_R = \sqrt{\frac{R}{t}}\) ): \[ \begin{align} \varepsilon&=\sigma_{R_B}^2 + \sigma_{R_M}^2=\frac{R_B}{t_B}+\frac{R_M}{t_M} = \frac{R_B}{t_B}+\frac{R_M}{T-t_B} \left( \equiv \frac{R_B}{T-t_M}+\frac{R_M}{t_M} \right) \\ \\ \frac{d \varepsilon}{d t_B} &= -\frac{R_B}{t_B^2} +\frac{R_M}{\left(T-t_B\right)^2} = -\frac{R_B}{t_B^2} + \frac{R_M}{t_M^2} = 0 \Rightarrow \frac{R_B}{t_B^2} = \frac{R_M}{t_M^2} \\ \\ \Rightarrow r = \frac{t_M}{t_b}&=\sqrt{\frac{R_M}{R_B}} \end{align} \] Note, the derivative could have been taken with respect to \(t_M\) using \(\varepsilon = \frac{R_B}{T-t_M}+\frac{R_M}{t_M} \) to achieve the same result.

The above ratio can be used to find the fraction of a total fixed time, \(T\), that should be apportioned to monitoring the source plus background, \(f_M\), and the fraction that should be apportioned to measuring only the background, \(f_B\). Noting that \(r = \frac{t_M}{t_b} \Leftrightarrow t_M=rt_B\) \[ \begin{align} f_B&=\frac{t_B}{T}=\frac{t_B}{t_B+t_M}=\frac{t_B}{t_B+rt_B}=\frac{1}{1+r} \Rightarrow t_B=\frac{T}{1+r} \\ \\ f_M&=\frac{t_M}{T}=\frac{t_M}{t_B+t_M} = \frac{t_M}{\frac{t_M}{r}+t_M} = \frac{r}{1+r} \Rightarrow t_M=\frac{rT}{1+r} \end{align} \]

Using the values in Example 2 the optimal timing ratio for measuring the source with background to measuring just the background is \[ r = \frac{t_M}{t_b} = \sqrt{\frac{R_M}{R_B}} = \sqrt{\frac{50}{30}} = 1.291 \] giving \[ f_M = \frac{r}{1+r} = \frac{1.291}{2.291} = 0.56 \] \[ f_B = \frac{1}{1+r} = \frac{1}{2.291} = 0.44 \] i.e. you should spend 56% of the monitoring time on the source with background and 44% of the monitoring time obtaining a background value.

To explore this using the simulator change the operation mode in the lower left from "Independent times" to "Fixed total time". You can now give a total monitoring time and change the fraction of theat time that takes measurements of the source (plus background). Moving the fraction slider will trigger a resampling and you can see how the constructed rate error changes as more or less time is spent monitoring the source.

Note, in the above the constructed source error may jump around due to the random nature of the problem. However, if you collect enough counts you will see the minimum error is at a source sampling fraction time of 0.56. You can do this in the simulator by increasing the total smapling time (which reduces the errors). For the maximum total sampling time of 1440 minutes (24 hours) a simple swipe of the fraction slider allows you to see where the optimum should be.

Activity Change the actual source rate to 150 cpm and set the total sampling time to 1440 minutes. Use the fraction slider to estimate the optimum split and check your answer against the value calculated using the above result of \(f_M = \frac{r}{1+r}\) where \(r=\sqrt{\frac{R_M}{R_B}}\) - remember \(R_M\) is the combined rate of the source and background! You can take the actual values (assume perfect knowledge), or the sampled values to simulate how you might estimate a split based on a set of initial survey measurements made over a shorter time. Note, the errors in these shorter measurements might mean you don't get the theoretically perfect split - the question is how much does that matter?

What if you don't get the timing split right?

The figure below shows how the split of times between monitoring the source with background and monitoring the just the background affects the total source-only error for the values used in Example 2. As can be seen for the count rates chosen the total error is fairly flat over a wide range of time splits.

For the case where the \(R_M= 5000\) rather than 50, the following error sensitivity is seen. In this case you probably want to be nearer the optimal split.

Obtaining a specified error level - answered!

Under the assumption that the monitoring is going to be carried out using the optimal time split between monitoring the source with background and monitoring the just the background the question of what length of monitoring time is required to obtain a given accuracy when two measurements are required. As before the percentage error on \(R_S\) is \[\% \sigma_{R_S} = \frac{100 \sigma_{R_S}}{R_S} \Rightarrow \sigma_{R_S} = \frac{\% \sigma_{R_S} R_S}{100} \] Squaring both sides and using \(\sigma_{R_S} = \sqrt{\sigma_{R_B}^2 + \sigma_{R_M}^2}\) and \(\sigma_R = \sqrt{\frac{R}{t}}\) gives \[ {\left( \frac{\% \sigma_{R_S} R_S}{100} \right)}^2 = \sigma_{R_B}^2 + \sigma_{R_M}^2 = \frac{R_B}{t_B}+\frac{R_M}{t_M} = \frac{t_M R_B + t_B R_M}{t_B t_M} \] Noting \(t_M=r t_B\) and \(R_S = R_M - R_B\) \[ \begin{align} {\left( \frac{\% \sigma_{R_S} \left( R_M - R_B \right)}{100} \right)}^2 &= \frac{t_M R_B + t_B R_M}{t_B t_M} \\ \\ &=\frac{t_B R_M + r t_B R_B}{r t_B^2} \\ \\ &=\frac{R_M+r R_B}{r t_B} \end{align} \] Rearranging gives \[ t_B=\frac{{100}^2 \left( R_M + r R_B\right)}{r {\left( \% \sigma_{R_S} \left( R_M - R_B \right) \right)}^2} \] And \[ t_M = rt_B = \frac{{100}^2 \left( R_M + r R_B\right)}{{\left( \% \sigma_{R_S} \left( R_M - R_B \right) \right)}^2} \]

For the values in Example 2, \( R_B=30, R_M=50\) and \( r = \frac{t_M}{t_b} = \sqrt{\frac{R_M}{R_B}} = \sqrt{\frac{50}{30}} = 1.291\), the following times are obtained \[ \begin{align} t_B&=\frac{{100}^2 \left( 50 + 1.291 \times 30\right)}{1.291 \times {\left( 1 \times \left( 50 - 30 \right) \right)}^2} \approx 1718 \space minutes \\ \\ t_M&=\frac{{100}^2 \left( 50 + 1.291 \times 30\right)}{{\left( 1 \times \left( 50 - 30 \right) \right)}^2} \approx 2218 \space minutes \\ \\ \end{align} \] The total monitoring time is \(T=t_B+t_M \approx 3936\) minutes (65.6 hours or 2.73 days).

Recall that for a single measurement the time required to obtain a given percentage error is \(t= \frac{1}{R} {\left( \frac{100}{\% \sigma_R} \right)}^2\). Individually for \(R_B = 30\) and \(R_M = 50\) the time required to achieve a 1% error is \[ \begin{align} t_B^{ind.}&=\frac{1}{30}{\left( \frac{100}{1} \right)}^2 = 333 \space minutes \\ \\ t_M^{ind.}&=\frac{1}{50}{\left( \frac{100}{1} \right)}^2 = 200 \space minutes \end{align} \] Note, these times are an order of magnitude shorter than the time required to achieve a 1% error when two independent measurements are made and used together to estimate the source only rate.

Is the source really there?

When the activity of the source is large compared with the background (or even comparable) it is likely that monitoring will confirm the presence of a source. However, when the source activity is small, or there is no source can I (a) detect it or (b) definitively say it's not there?

The answer to the above is related to testing errors:

• Obtaining a false positive, i.e. saying there is a source present when in fact there isn't. This is a Type I error.
• Obtaining a false negative, i.e. saying there is no source present when in fact there is. This is a Type II error.

Type I error - a false positive

In the case of a false positive, the expected value of \(R_S = R_M - R_B\) is zero as there is no source. However, as the constructed value of \(R_S\) has been made with two independent measurements it has a combined standard deviation of \(\sigma_{R_S} = \sqrt{\sigma_{R_B}^2 + \sigma_{R_M}^2}\). An example will help with the discussion.

Assume the following results were obtained when measuring \(R_B\) and \(R_M\): \[ \begin{align} R_B &= 8.100 \pm 0.900 \\ \\ R_M &= 11.200 \pm 1.058 \end{align} \] These give a constructed source rate of \[ R_S = 3.100 \pm 1.389 \] Approximating the Poisson distributions with normal distributions (careful!), the value of \(R_S\) is distributed about zero with a standard deviation of 1.389.

Notes:

• It is centred on zero as we know there is no source so with perfect measurement \(R_B = R_M \Rightarrow R_S=0\).
• The constructed value of \(R_S\) can be less than zero as \(R_B\) and \(R_S\) are independly taken from the same distribution and there will be occasions where \(R_B \gt R_M\).
• The standard deviation from \(R_S\) is used as this is our combined measurement standard deviation.

The screen shot below shows an example of the second bullet point above in the simulator.

To determine whether the value \(R_S=3.1\) represents an actual source or is within measurement uncertainty the one-sided right hand tail is examined. The right-tailed critical values for various significance levels are

If you were to say anything that has less than a 10% or 5% chance of occurring at random is regarded as a true signal you would conclude that as \(R_S \gt1.778\) or \(R_S \gt 2.285\) then there is a source present with the given rate. This is a false positive reading. However, if you say that anything that has a less than 1% chance of occurring at random is regarded as a true signal you would conclude that there is no source present. This latter conclusion may be problematic for cases where there is actually a source present and it has been rejected as this leads to the Type II error - a false negative.

Type II error - a false negative

Suppose now there is a source with a count rate of 3.1 cpm. The actual value detected in the absence of any background will follow a distribution with a standard deviation depending on the length of the monitoring time. Superimposing this distribution over the distribution of an assumed absence of a source gives

The dotted lines in this figure show, from left to right, the cut off measured values of \(R_S\) to avoid a false positive at the 10%, 5% and 1% levels. However, it is seen that these lines also significantly overlap the distribution of there actually being a source present. In this example, avoiding a Type I error at the 10% level (the area under the blue curve to the right of the 10% dotted line ) means there is a 17% chance that you have missed an actual source being present (the area under the orange curve to the left of the 10% dotted line) - i.e. a 17% chance of making a type II error. Things get worse when you further decrease the chance of making a Type I error - see below.

The Type II significance is the probability of making a Type II (false negative) error. The right hand column is one minus the Type II significance and it gives the probability that you can correctly detect the presence of the source.

Controlling Type I and Type II errors

The above has highlighted Type I (false positive) and Type II (false negative) errors. In relation to detection:

• The controlling of Type I errors is related to the minimum significant activity that can be detected.
• The controlling of Type II errors is related to the minimum true activity that can be detected.

As Type I and Type II errors are dependent on the standard deviation of a measurement which, in a single measurement case, is determined by the monitoring time, the way to control the errors is to ensure a sufficient monitoring time has been taken. The figure below shows the Type I/II figure above but with distributions that have smaller standard deviations (as a result of increased monitoring time).

In the above, allowing for a 1% chance of a Type I error leads to only a 3.6% chance of a Type II error.