# Comparison of Two Approaches for Detection and Estimation of Radioactive Sources

Authors: Eugene Yee; Ajith Gunatilaka; Branko Ristic
Journal: ISRN Applied Mathematics
 (2011)

---

## Abstract

This paper describes and compares two approaches for the problem of determining the number of radioactive point sources that potentially exist in a designated area and estimating the parameters of these sources (their locations and strengths) using a small number of noisy radiological measurements provided by a radiation sensor. Both approaches use the Bayesian inferential methodology but sample the posterior distribution differently: one approach uses importance sampling with progressive correction and the other a reversible-jump Markov chain Monte Carlo sampling. The two approaches also use different measurement models for the radiation data. The first approach assumes a perfect knowledge of the data model and the average background radiation level, whereas the second approach quantifies explicitly the uncertainties in the model specification and in the average background radiation level. The performances of the two approaches are compared using experimental data acquired during a recent radiological field trial.


---

## Body

## 1. Introduction

There is growing concern in recent years regarding the risk of the smuggling and illicit trafficking of stolen radiological material. This state of affairs has heightened the spectre of radiological terrorism involving the use of an improvised radiological dispersal device (dirty bomb) to disperse radiological material over a large area using the force of conventional explosives (see, Panofsky [1]). In this context, the detection, localisation, and characterization of sources of radioactive material using radiation sensors (or, networks of such sensors) have important implications for safety and public security.Various researchers have focussed on the problem of recovering the physical location and strength of radioactive sources using measurements obtained from a number of radiological detectors. Howse et al. [2] applied a nonlinear, recursive least-squares algorithm for tracking the positions of a moving radioactive point source. Stephens Jr. and Peurrung [3] conducted a cost-benefit analysis involving a variety of (conflicting) operational requirements that need to be considered in the design and deployment of a network of radiation sensors for use in the detection of a moving radioactive source. Nemzek et al. [4] and Brennan et al. [5] investigated the application of a Bayesian methodology for the detection and estimation of mobile radioactive sources, with a specific focus on examination of the signal-to-noise ratio envelope (and the concomitant detection limits) expected from utilization of a distributed sensor network for detection of a radiological source in a moving vehicle. Gunatilaka et al. [6] proposed the application of a maximum likelihood estimation algorithm, extended Kalman filter, and unscented Kalman filter for the localisation of a single static radiological point source. Morelande et al. [7, 8] applied a Bayesian methodology realized using a particle filtering technique with progressive correction to localize multiple static radiological point sources. Finally, Mendis et al. [9] applied binary and continuous genetic algorithms to solve a maximum likelihood estimation problem for recovery of the location and strengths of an a priori known number of static radiological point sources.In this paper, we describe two approaches for radiological source localisation for the case when the number of sources is unknowna priori and compare their performance by applying them to an experimental data set collected in a field trial in Puckapunyal Military Area (Victoria, Australia). Both approaches use the Bayesian inferential methodology but sample the posterior distribution using different algorithms. In the first approach, the sampling algorithm uses importance sampling with progressive correction (PC), while the second approach uses the Markov chain Monte Carlo (MCMC) algorithm for the sampling. While the first algorithm uses the minimum description length (MDL) to estimate the unknown number of sources, the second algorithm uses the reversible-jump MCMC approach for this purpose. Finally, the measurement model for the radiation data differs for the two approaches in the following way: the first approach uses a form of the posterior distribution that implicitly assumes that the average background radiation count rate is exactly known and that the error in the model for the radiation data is negligible, whereas the second approach utilizes a form of the posterior distribution that relaxes these assumptions.The experimental data for evaluating the two approaches were acquired using a Geiger-Müller detector in the presence of multiple radiological point sources of unequal strengths. Background radiation data in the absence of these radiation sources were also acquired to provide an estimate of the average background radiation level.The organisation of the paper is as follows. The radiological source localisation problem is formulated in Section2. Section 3 describes the two approaches we use for radiological source estimation. Section 4 describes the radiation field trial and experimental data sets. Some implementation aspects of the two algorithms are discussed in Section 5. Section 6 presents the results of the application of the two approaches for source reconstruction to the experimental data. Conclusions are drawn in Section 7.


## 2. Problem Formulation

The problem of the detection and estimation of the characteristics (e.g., location, activity) of an unknown number of static radiological point sources from a finite number of “noisy” measurements (data) obtained from radiation sensors is addressed using the Bayesian formalism. This formalism involves the application of Bayes' rule:(2.1)p(Θ∣z,I)∝p(Θ∣I)p(z∣Θ,I),
where Θ are the parameters that we are trying to determine and z is the data. In  (2.1), I is the background (contextual) information available in the problem, “∣” denotes “conditional upon”, and the various factors that appear in this relation have the following interpretation: p(Θ∣I) is the prior probability for the parameters Θ, p(z∣Θ,I) is the likelihood function which specifies the probability that we observe data z when Θ is known exactly, and p(Θ∣z,I) is the posterior probability for the parameters Θ in light of the new information introduced through the acquired data z. In essence, Bayes' rule informs us how to update our knowledge about the parameters Θ encoded in p(Θ∣I) following the acquisition of data z. To apply the Bayesian formalism to the radiological source reconstruction problem, we need to assign appropriate functional forms for the prior probability and the likelihood function, which in turn determines the posterior probability.


### 2.1. Assignment of Prior Probability

Here, we assume that a numberr of point sources of gamma radiation are present in a flat open area that is free of any obstructions. The number of sources, r, is unknown. If  r>0, the ith point source of gamma radiation (where i=1,2,…,r) is parameterised by  (i)
its locationxs,i≡(xs,i,ys,i) in a Cartesian coordinate system with (xs,i,ys,i)∈𝒜, where 𝒜⊂ℝ2 is some large region that is assumed to contain all the sources,(ii)
its intensity rate (source strength),Qs,i.The parameter vector of sourcei is thus given by(2.2)θi=(xs,iys,iQs,i)⊺,
where ⊺ denotes the matrix transpose. The source parameter vectors are collected into a stacked vector θ(r)=(θ1⊺⋯θr⊺)⊺, and let Θ≡(rθ(r)⊺)⊺.If we assume the logical independence of the various source parameters (or equivalently, the components ofΘ), the prior probability for the parameters factorizes as follows:(2.3)p(Θ∣I)=p(r∣I)∏i=1rp(Qs,i∣I)p(xs,i∣I).
Now, we need to assign explicit functional forms for each of the component prior probabilities in  (2.3).The prior distribution ofr,  p(r∣I), is assigned as a binomial distribution:(2.4)p(r∣I)=(rmax⁡-rmin⁡)!(r-rmin⁡)!(rmax⁡-r)!p*(r-rmin⁡)(1-p*)rmax⁡-r,
(r=rmin,rmin+1,…,rmax). In  (2.4), p*∈(0,1) is the binomial rate, and, rmin and rmax are, respectively, the minimum and maximum number of sources.Two alternative functional forms are used for the prior distribution ofQs,i (i=1,2,…,r). The first of the two alternatives assigns the prior for the source strength Qs,i to be a gamma distribution:(2.5)p(Qs,i∣I)=Qs,iκ-1exp⁡⁡(-Qs,i/ψ)Γ(κ)ψκ,i=1,2,…,r,
where κ and ψ are the shape and scale parameters, respectively, and Γ(x) is the gamma function. For the second of the two alternatives, the prior distribution for Qs,i is assigned a Bernoulli-uniform mixture model:(2.6)p(Qs,i∣I)=(1-γ)δ(Qs,i)+γI(0,Qmax⁡)(Qs,i)Qmax⁡,i=1,2,…,r,
where γ is the probability that the source is active (Pr{Qs,i>0}=γ), δ(x) is the Dirac delta function, and Qmax is an a priori upper bound on the expected source intensity rate. In  (2.6), 𝕀A(x) denotes the indicator (characteristic) function for set A, with 𝕀A(x)=1 if x∈A and 𝕀A(x)=0 if x∉A.Finally, the prior distribution for the source locationxs,i is assigned to be uniformly distributed over the region 𝒜⊂ℝ2 that is assumed a priori to contain the source:(2.7)p(xs,i∣I)=IA(xs,i)meas⁡(A),i=1,2,…,r,
where meas(𝒜) is the area of the region 𝒜.


### 2.2. Assignment of Likelihood Function

To prescribe a functional form for the likelihood function, we need to relate the hypotheses of interest about the unknown source(s)1 as encoded in Θ to the available radiation data z measured by detectors emplaced in the vicinity of the source(s). Towards this objective, we need to formulate a measurement model for the radiation data z.The radiation counts from nuclear decay obey Poisson statistics (see, Tsoulfanidis [10]), implying that the probability that a gamma radiation detector registers z∈ℕ∪{0} counts (ℕ∪{0} being the set of natural numbers including zero) in an exposure time of τ seconds, from a source that emits on average μ counts per second, is(2.8)P(z;λ)=λzz!exp⁡⁡(-λ),
where λ=μτ is the parameter (both the mean and the variance) of the Poisson distribution.The measurements of the radiation field are made using a Geiger-Müller (GM) counter. Letzj∈ℕ∪{0} (j=1,…,m) be a measured count from the GM counter, taken at (receptor) location (xl,j,yl,j)∈ℝ2. The following assumptions are made: (1) the GM counter has a uniform directional response; (2) the radiation sources are point sources; (3) there is negligible attenuation of gamma radiation due to air; (4) the radiation measurements are independently distributed, and, (5) the exposure time (or the sampling interval) τ is constant for all measurements.


#### 2.2.1. Measurement Model: Known Background Rate

The likelihood function is simply the joint density of the measurement vectorz≡(z1⋯zm)⊺ conditional on the parameter vector θ and the knowledge that r sources are present. With this identification, the likelihood function can then be written as a product of Poisson distributions (Martin and Harbison [11]; Tsoulfanidis [10]) as follows:(2.9)p(z∣μb,Θ,I)≡p(z∣μb,r,θ,I)≡l(z∣μb,r,θ)=∏j=1mP(zj;λj(Θ)),
where λj(Θ) is the mean radiation count for the jth sensor location:(2.10)λj(Θ)=(μb+∑i=1rQs,idji2)τ,
with(2.11)dji=((xl,j-xs,i)2+(yl,j-ys,i)2)1/2
being the Euclidean distance between the ith source and the jth sensor.The constantμb in  (2.10) is the mean background count rate (namely, the average count rate due to the background radiation only which includes the contribution from cosmic and terrestrial radiation). In the measurement model (and concomitant form of the likelihood function given by  (2.9) and (2.10)), the mean background count rate μb and the number of radiological point sources r are assumed to be known a priori. In consequence, μb and r have been added as quantities to the right of the vertical bar in the likelihood function of  (2.9) to indicate explicitly that these quantities are known exactly.The mean signal count rateμsj at the jth sensor due to the radioactive point sources only is modelled in  (2.10) as(2.12)μ̃sj=∑i=1rQs,i(xl,j-xs,i)2+(yl,j-ys,i)2.
In  (2.10), it is implicitly assumed that the model for the mean signal count rate in  (2.12) is perfect (no model error) so μsj=μ̃sj.


#### 2.2.2. Measurement Model: Unknown Background Rate

The key assumptions in the formulation of the measurement model of  (2.9) and (2.10) are that (1) the mean background count rate μb is known exactly, and (2) the model error associated with the model for the mean signal count rate μs of  (2.12) is negligible. Now, we will relax these critical assumptions and approach the (perhaps) more realistic case when we must incorporate both the model error and the effects of an uncertain mean background count rate in the measurement model.Towards this objective, the following model is assumed for the (true though unknown) mean signal count rateμsj (j=1,2,…,m). The mean signal count rate can be described by the radiation data z and by the quantitative model (cf.  (2.12)), so
(2.13)μsj=μ̂sj+ej(1)=μ̃sj+ej(2),j=1,2,…,m,
where μ̂sj is an estimate for the mean signal count rate at the jth sensor obtained from the data z (more specifically, zj), μ̃sj is the mean signal count rate at the jth sensor obtained from the model given in (2.12), ej(1) is the uncertainty in the estimate μ̂sj of the mean signal count rate, and ej(2) represents the model error incurred by using (2.12) to predict the mean signal count rate. The model errors consist of a number of contributions which include (1) the effects of absorption of the gamma radiation by the air and (2) the effects of backscattering from the ground surface which may (possibly) influence the inverse-square law dependence in the model of (2.12).The measurement model of  (2.13) can be rearranged to give(2.14)μ̂sj-μ̃sj=ej(2)-ej(1)≡ϵj,j=1,2,…,m,
where ϵj is the composite error. If the expectation value 〈ϵj〉 of the composite error is zero and the variance 〈ϵj2〉 is assumed to be given by σj2, then application of the principle of maximum entropy (see, Jaynes [12]) provides the following explicit form for the likelihood function (assuming that at least one radiological source is present, so r>0):(2.15)p(z∣Θ,I)=1∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}.
It is noted that the Gaussian likelihood (which is effectively the probability distribution of the “noise” ϵj) obtained from application of the maximum entropy principle is simply the most conservative noise distribution that is consistent with the given facts (namely, the first two moments of the “noise” are known). In the Gaussian likelihood function of (2.15), the variance of the composite error ϵj is obtained as(2.16)σj2=σ1,j2+σ2,j2,j=1,2,…,m,
where σ1,j2 is the variance of the uncertainty ej(1) in the estimate μ̂sj of the mean signal count rate (obtained from the data z) and σ2,j2 is the variance of the model error ej(2). In this paper, the model error standard deviation (or square root of the variance) is specified as σ2,j=(μ̃sj)1/2.To complete the specification for the likelihood function in  (2.15), we need to provide the estimate μ̂sj along with the variance σ1,j2 in this estimate. To this purpose, we assume that we have available a measurement of background radiation (obtained in the absence of any point radiation sources) given to us as zb counts obtained over an exposure time τb. Now, suppose that a radiation sensor at a fixed receptor location obtains a measurement of z counts in an exposure time τ. We will now show how to use this information to obtain an estimate for the mean signal count rate μ̂s, as well as to characterize the uncertainty in this estimate in terms of the standard deviation σ1.To accomplish this, we utilize a Bayesian solution for the extraction of weak signals in strong background interference described by Loredo [13]. We will briefly outline the application of this methodology to the problem in this paper. To begin, we can apply Bayes' rule to compute the joint probability for the mean signal count rate μs and the mean background count rate μb conditioned on the measurement z. This gives(2.17)p(μs,μb∣z,I)∝p(μs∣μb,I)p(μb∣I)p(z∣μs,μb,I).
It is stressed that the contextual information I in (2.17) includes the contextual information Ib concerning the measurement of the background radiation and the contextual information Is concerning the measurement of the radiation field which may include the presence of one or more radiation sources. Hence, the contextual information I that is available in the problem can be decomposed as I=Ib,Is.We will assign a least informative prior probability forμs (the so-called Jeffreys prior; see Jaynes [12]) which has the following form (conditioned on a known mean background rate μb):(2.18)p(μs∣μb,I)=1μs+μb.
The prior probability for μb is chosen to be the following informative prior (using the information that is available, namely, that a measurement of the background radiation yielded zb counts in an exposure time τb):(2.19)p(μb∣I)≡p(μb∣zb,Ib,Is)=p(μb∣zb,Ib)∝p(μb∣Ib)p(zb∣μb,Ib),
since Is does not have any bearing on the determination of μb. The prior p(μb∣Ib) is assigned the nonsinformative Jeffreys prior so p(μb∣Ib)=1/μb and the likelihood function p(zb∣μb,Ib)=(μbτb)zbexp(-μbτb)/zb!, so(2.20)p(μb∣I)=p(μb∣zb,Ib)∝τb(μbτb)zb-1exp⁡⁡(-μbτb).The likelihood function in (2.17), associated with the measurement of z counts in exposure interval τ coming from a radiation field with mean count rate (μs+μb), is described by the following Poisson distribution:(2.21)p(z∣μs,μb,I)=τz(μs+μb)zz!exp⁡⁡(-(μs+μb)τ).
Now, if we substitute (2.18), (2.20), and (2.21) in  (2.17), the joint probability for the mean signal and background count rates is given by(2.22)p(μs,μb∣z,I)∝(μs+μb)z-1μbzb-1exp⁡⁡(-μsτ)exp⁡⁡(-μb(τb+τ)).
The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the (unknown) mean background count rate:(2.23)p(μs∣z,I)=∫0∞p(μs,μb∣z,I)dμb,
which on using  (2.22) and performing the integration finally gives (on normalization of the probability density function)(2.24)p(μs∣z,I)=∑i=1zΦiτ(μsτ)i-1(i-1)!exp⁡⁡(-μsτ),
where(2.25)Φi=(1+τb/τ)i(z+zb-i-1)!/(z-i)!∑j=1z(1+τb/τ)j(z+zb-j-1)!/(z-j)!.
Note that  (2.24) and (2.25) provide an explicit expression for the posterior probability of the mean signal count rate μs that is independent of the (unknown) mean background count rate μb. More specifically, the precise value of μb is not known, so it is treated as a nuisance parameter here and removed by integration (marginalization).Given the posterior probability ofμs summarized by (2.24) and (2.25), a best estimate for the mean signal count rate can be obtained as the posterior mean:(2.26)μ̂s≡〈μs〉=∫0∞μsp(μs∣z,I)dμs=1τ∑i=1ziΦi.
The second moment of μs about zero can be evaluated as(2.27)〈μs2〉=∫0∞μs2p(μs∣z,I)dμs=1τ2∑i=1zi(i+1)Φi,
from which the uncertainty σ1 in the estimate for μs given in  (2.26) can be determined as(2.28)σ1=(〈μs2〉-(μ̂s)2)1/2.
Equations (2.26), (2.27) and (2.28) when used in conjunction with (2.16) and (2.15) fully determine the form of the likelihood function (assuming implicitly that r>0). This functional form for the likelihood function accounts explicitly for the model error and the sampling error (the latter of which is compounded by an uncertain mean background count rate).For the caser=0, the likelihood function needs to be determined as follows (again assuming that our state of knowledge consists of a background radiation measurement yielding zb counts in exposure time τb):(2.29)p(z∣zb,I)=∏j=1m∫0∞P(zj∣μbτ)p(μb∣zb,Ib)dμb,
which on insertion of the normalized form of  (2.20) for p(μb∣zb,Ib) and evaluation of the integral gives(2.30)p(z∣zb,I)=∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!.
In (2.29) and (2.30), the radiation measurements zj are assumed to be made using the same exposure time τ (without any loss in generality).Finally, the likelihood function of (2.15) valid for r>0 and the likelihood function of (2.30) valid for r=0 can be combined to give the following general likelihood function for the measurement model assumed in this section:(2.31)p(z∣zb,Θ,I)=I{r>0}(r)∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}+I{r=0}(r)∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!,
where the assumption of the availability of a background radiation count zb (obtained over an exposure time τb) has been included in the conditioning of the likelihood function (and of the posterior distribution based upon this likelihood function).


## 2.1. Assignment of Prior Probability

Here, we assume that a numberr of point sources of gamma radiation are present in a flat open area that is free of any obstructions. The number of sources, r, is unknown. If  r>0, the ith point source of gamma radiation (where i=1,2,…,r) is parameterised by  (i)
its locationxs,i≡(xs,i,ys,i) in a Cartesian coordinate system with (xs,i,ys,i)∈𝒜, where 𝒜⊂ℝ2 is some large region that is assumed to contain all the sources,(ii)
its intensity rate (source strength),Qs,i.The parameter vector of sourcei is thus given by(2.2)θi=(xs,iys,iQs,i)⊺,
where ⊺ denotes the matrix transpose. The source parameter vectors are collected into a stacked vector θ(r)=(θ1⊺⋯θr⊺)⊺, and let Θ≡(rθ(r)⊺)⊺.If we assume the logical independence of the various source parameters (or equivalently, the components ofΘ), the prior probability for the parameters factorizes as follows:(2.3)p(Θ∣I)=p(r∣I)∏i=1rp(Qs,i∣I)p(xs,i∣I).
Now, we need to assign explicit functional forms for each of the component prior probabilities in  (2.3).The prior distribution ofr,  p(r∣I), is assigned as a binomial distribution:(2.4)p(r∣I)=(rmax⁡-rmin⁡)!(r-rmin⁡)!(rmax⁡-r)!p*(r-rmin⁡)(1-p*)rmax⁡-r,
(r=rmin,rmin+1,…,rmax). In  (2.4), p*∈(0,1) is the binomial rate, and, rmin and rmax are, respectively, the minimum and maximum number of sources.Two alternative functional forms are used for the prior distribution ofQs,i (i=1,2,…,r). The first of the two alternatives assigns the prior for the source strength Qs,i to be a gamma distribution:(2.5)p(Qs,i∣I)=Qs,iκ-1exp⁡⁡(-Qs,i/ψ)Γ(κ)ψκ,i=1,2,…,r,
where κ and ψ are the shape and scale parameters, respectively, and Γ(x) is the gamma function. For the second of the two alternatives, the prior distribution for Qs,i is assigned a Bernoulli-uniform mixture model:(2.6)p(Qs,i∣I)=(1-γ)δ(Qs,i)+γI(0,Qmax⁡)(Qs,i)Qmax⁡,i=1,2,…,r,
where γ is the probability that the source is active (Pr{Qs,i>0}=γ), δ(x) is the Dirac delta function, and Qmax is an a priori upper bound on the expected source intensity rate. In  (2.6), 𝕀A(x) denotes the indicator (characteristic) function for set A, with 𝕀A(x)=1 if x∈A and 𝕀A(x)=0 if x∉A.Finally, the prior distribution for the source locationxs,i is assigned to be uniformly distributed over the region 𝒜⊂ℝ2 that is assumed a priori to contain the source:(2.7)p(xs,i∣I)=IA(xs,i)meas⁡(A),i=1,2,…,r,
where meas(𝒜) is the area of the region 𝒜.


## 2.2. Assignment of Likelihood Function

To prescribe a functional form for the likelihood function, we need to relate the hypotheses of interest about the unknown source(s)1 as encoded in Θ to the available radiation data z measured by detectors emplaced in the vicinity of the source(s). Towards this objective, we need to formulate a measurement model for the radiation data z.The radiation counts from nuclear decay obey Poisson statistics (see, Tsoulfanidis [10]), implying that the probability that a gamma radiation detector registers z∈ℕ∪{0} counts (ℕ∪{0} being the set of natural numbers including zero) in an exposure time of τ seconds, from a source that emits on average μ counts per second, is(2.8)P(z;λ)=λzz!exp⁡⁡(-λ),
where λ=μτ is the parameter (both the mean and the variance) of the Poisson distribution.The measurements of the radiation field are made using a Geiger-Müller (GM) counter. Letzj∈ℕ∪{0} (j=1,…,m) be a measured count from the GM counter, taken at (receptor) location (xl,j,yl,j)∈ℝ2. The following assumptions are made: (1) the GM counter has a uniform directional response; (2) the radiation sources are point sources; (3) there is negligible attenuation of gamma radiation due to air; (4) the radiation measurements are independently distributed, and, (5) the exposure time (or the sampling interval) τ is constant for all measurements.


### 2.2.1. Measurement Model: Known Background Rate

The likelihood function is simply the joint density of the measurement vectorz≡(z1⋯zm)⊺ conditional on the parameter vector θ and the knowledge that r sources are present. With this identification, the likelihood function can then be written as a product of Poisson distributions (Martin and Harbison [11]; Tsoulfanidis [10]) as follows:(2.9)p(z∣μb,Θ,I)≡p(z∣μb,r,θ,I)≡l(z∣μb,r,θ)=∏j=1mP(zj;λj(Θ)),
where λj(Θ) is the mean radiation count for the jth sensor location:(2.10)λj(Θ)=(μb+∑i=1rQs,idji2)τ,
with(2.11)dji=((xl,j-xs,i)2+(yl,j-ys,i)2)1/2
being the Euclidean distance between the ith source and the jth sensor.The constantμb in  (2.10) is the mean background count rate (namely, the average count rate due to the background radiation only which includes the contribution from cosmic and terrestrial radiation). In the measurement model (and concomitant form of the likelihood function given by  (2.9) and (2.10)), the mean background count rate μb and the number of radiological point sources r are assumed to be known a priori. In consequence, μb and r have been added as quantities to the right of the vertical bar in the likelihood function of  (2.9) to indicate explicitly that these quantities are known exactly.The mean signal count rateμsj at the jth sensor due to the radioactive point sources only is modelled in  (2.10) as(2.12)μ̃sj=∑i=1rQs,i(xl,j-xs,i)2+(yl,j-ys,i)2.
In  (2.10), it is implicitly assumed that the model for the mean signal count rate in  (2.12) is perfect (no model error) so μsj=μ̃sj.


### 2.2.2. Measurement Model: Unknown Background Rate

The key assumptions in the formulation of the measurement model of  (2.9) and (2.10) are that (1) the mean background count rate μb is known exactly, and (2) the model error associated with the model for the mean signal count rate μs of  (2.12) is negligible. Now, we will relax these critical assumptions and approach the (perhaps) more realistic case when we must incorporate both the model error and the effects of an uncertain mean background count rate in the measurement model.Towards this objective, the following model is assumed for the (true though unknown) mean signal count rateμsj (j=1,2,…,m). The mean signal count rate can be described by the radiation data z and by the quantitative model (cf.  (2.12)), so
(2.13)μsj=μ̂sj+ej(1)=μ̃sj+ej(2),j=1,2,…,m,
where μ̂sj is an estimate for the mean signal count rate at the jth sensor obtained from the data z (more specifically, zj), μ̃sj is the mean signal count rate at the jth sensor obtained from the model given in (2.12), ej(1) is the uncertainty in the estimate μ̂sj of the mean signal count rate, and ej(2) represents the model error incurred by using (2.12) to predict the mean signal count rate. The model errors consist of a number of contributions which include (1) the effects of absorption of the gamma radiation by the air and (2) the effects of backscattering from the ground surface which may (possibly) influence the inverse-square law dependence in the model of (2.12).The measurement model of  (2.13) can be rearranged to give(2.14)μ̂sj-μ̃sj=ej(2)-ej(1)≡ϵj,j=1,2,…,m,
where ϵj is the composite error. If the expectation value 〈ϵj〉 of the composite error is zero and the variance 〈ϵj2〉 is assumed to be given by σj2, then application of the principle of maximum entropy (see, Jaynes [12]) provides the following explicit form for the likelihood function (assuming that at least one radiological source is present, so r>0):(2.15)p(z∣Θ,I)=1∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}.
It is noted that the Gaussian likelihood (which is effectively the probability distribution of the “noise” ϵj) obtained from application of the maximum entropy principle is simply the most conservative noise distribution that is consistent with the given facts (namely, the first two moments of the “noise” are known). In the Gaussian likelihood function of (2.15), the variance of the composite error ϵj is obtained as(2.16)σj2=σ1,j2+σ2,j2,j=1,2,…,m,
where σ1,j2 is the variance of the uncertainty ej(1) in the estimate μ̂sj of the mean signal count rate (obtained from the data z) and σ2,j2 is the variance of the model error ej(2). In this paper, the model error standard deviation (or square root of the variance) is specified as σ2,j=(μ̃sj)1/2.To complete the specification for the likelihood function in  (2.15), we need to provide the estimate μ̂sj along with the variance σ1,j2 in this estimate. To this purpose, we assume that we have available a measurement of background radiation (obtained in the absence of any point radiation sources) given to us as zb counts obtained over an exposure time τb. Now, suppose that a radiation sensor at a fixed receptor location obtains a measurement of z counts in an exposure time τ. We will now show how to use this information to obtain an estimate for the mean signal count rate μ̂s, as well as to characterize the uncertainty in this estimate in terms of the standard deviation σ1.To accomplish this, we utilize a Bayesian solution for the extraction of weak signals in strong background interference described by Loredo [13]. We will briefly outline the application of this methodology to the problem in this paper. To begin, we can apply Bayes' rule to compute the joint probability for the mean signal count rate μs and the mean background count rate μb conditioned on the measurement z. This gives(2.17)p(μs,μb∣z,I)∝p(μs∣μb,I)p(μb∣I)p(z∣μs,μb,I).
It is stressed that the contextual information I in (2.17) includes the contextual information Ib concerning the measurement of the background radiation and the contextual information Is concerning the measurement of the radiation field which may include the presence of one or more radiation sources. Hence, the contextual information I that is available in the problem can be decomposed as I=Ib,Is.We will assign a least informative prior probability forμs (the so-called Jeffreys prior; see Jaynes [12]) which has the following form (conditioned on a known mean background rate μb):(2.18)p(μs∣μb,I)=1μs+μb.
The prior probability for μb is chosen to be the following informative prior (using the information that is available, namely, that a measurement of the background radiation yielded zb counts in an exposure time τb):(2.19)p(μb∣I)≡p(μb∣zb,Ib,Is)=p(μb∣zb,Ib)∝p(μb∣Ib)p(zb∣μb,Ib),
since Is does not have any bearing on the determination of μb. The prior p(μb∣Ib) is assigned the nonsinformative Jeffreys prior so p(μb∣Ib)=1/μb and the likelihood function p(zb∣μb,Ib)=(μbτb)zbexp(-μbτb)/zb!, so(2.20)p(μb∣I)=p(μb∣zb,Ib)∝τb(μbτb)zb-1exp⁡⁡(-μbτb).The likelihood function in (2.17), associated with the measurement of z counts in exposure interval τ coming from a radiation field with mean count rate (μs+μb), is described by the following Poisson distribution:(2.21)p(z∣μs,μb,I)=τz(μs+μb)zz!exp⁡⁡(-(μs+μb)τ).
Now, if we substitute (2.18), (2.20), and (2.21) in  (2.17), the joint probability for the mean signal and background count rates is given by(2.22)p(μs,μb∣z,I)∝(μs+μb)z-1μbzb-1exp⁡⁡(-μsτ)exp⁡⁡(-μb(τb+τ)).
The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the (unknown) mean background count rate:(2.23)p(μs∣z,I)=∫0∞p(μs,μb∣z,I)dμb,
which on using  (2.22) and performing the integration finally gives (on normalization of the probability density function)(2.24)p(μs∣z,I)=∑i=1zΦiτ(μsτ)i-1(i-1)!exp⁡⁡(-μsτ),
where(2.25)Φi=(1+τb/τ)i(z+zb-i-1)!/(z-i)!∑j=1z(1+τb/τ)j(z+zb-j-1)!/(z-j)!.
Note that  (2.24) and (2.25) provide an explicit expression for the posterior probability of the mean signal count rate μs that is independent of the (unknown) mean background count rate μb. More specifically, the precise value of μb is not known, so it is treated as a nuisance parameter here and removed by integration (marginalization).Given the posterior probability ofμs summarized by (2.24) and (2.25), a best estimate for the mean signal count rate can be obtained as the posterior mean:(2.26)μ̂s≡〈μs〉=∫0∞μsp(μs∣z,I)dμs=1τ∑i=1ziΦi.
The second moment of μs about zero can be evaluated as(2.27)〈μs2〉=∫0∞μs2p(μs∣z,I)dμs=1τ2∑i=1zi(i+1)Φi,
from which the uncertainty σ1 in the estimate for μs given in  (2.26) can be determined as(2.28)σ1=(〈μs2〉-(μ̂s)2)1/2.
Equations (2.26), (2.27) and (2.28) when used in conjunction with (2.16) and (2.15) fully determine the form of the likelihood function (assuming implicitly that r>0). This functional form for the likelihood function accounts explicitly for the model error and the sampling error (the latter of which is compounded by an uncertain mean background count rate).For the caser=0, the likelihood function needs to be determined as follows (again assuming that our state of knowledge consists of a background radiation measurement yielding zb counts in exposure time τb):(2.29)p(z∣zb,I)=∏j=1m∫0∞P(zj∣μbτ)p(μb∣zb,Ib)dμb,
which on insertion of the normalized form of  (2.20) for p(μb∣zb,Ib) and evaluation of the integral gives(2.30)p(z∣zb,I)=∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!.
In (2.29) and (2.30), the radiation measurements zj are assumed to be made using the same exposure time τ (without any loss in generality).Finally, the likelihood function of (2.15) valid for r>0 and the likelihood function of (2.30) valid for r=0 can be combined to give the following general likelihood function for the measurement model assumed in this section:(2.31)p(z∣zb,Θ,I)=I{r>0}(r)∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}+I{r=0}(r)∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!,
where the assumption of the availability of a background radiation count zb (obtained over an exposure time τb) has been included in the conditioning of the likelihood function (and of the posterior distribution based upon this likelihood function).


## 2.2.1. Measurement Model: Known Background Rate

The likelihood function is simply the joint density of the measurement vectorz≡(z1⋯zm)⊺ conditional on the parameter vector θ and the knowledge that r sources are present. With this identification, the likelihood function can then be written as a product of Poisson distributions (Martin and Harbison [11]; Tsoulfanidis [10]) as follows:(2.9)p(z∣μb,Θ,I)≡p(z∣μb,r,θ,I)≡l(z∣μb,r,θ)=∏j=1mP(zj;λj(Θ)),
where λj(Θ) is the mean radiation count for the jth sensor location:(2.10)λj(Θ)=(μb+∑i=1rQs,idji2)τ,
with(2.11)dji=((xl,j-xs,i)2+(yl,j-ys,i)2)1/2
being the Euclidean distance between the ith source and the jth sensor.The constantμb in  (2.10) is the mean background count rate (namely, the average count rate due to the background radiation only which includes the contribution from cosmic and terrestrial radiation). In the measurement model (and concomitant form of the likelihood function given by  (2.9) and (2.10)), the mean background count rate μb and the number of radiological point sources r are assumed to be known a priori. In consequence, μb and r have been added as quantities to the right of the vertical bar in the likelihood function of  (2.9) to indicate explicitly that these quantities are known exactly.The mean signal count rateμsj at the jth sensor due to the radioactive point sources only is modelled in  (2.10) as(2.12)μ̃sj=∑i=1rQs,i(xl,j-xs,i)2+(yl,j-ys,i)2.
In  (2.10), it is implicitly assumed that the model for the mean signal count rate in  (2.12) is perfect (no model error) so μsj=μ̃sj.


## 2.2.2. Measurement Model: Unknown Background Rate

The key assumptions in the formulation of the measurement model of  (2.9) and (2.10) are that (1) the mean background count rate μb is known exactly, and (2) the model error associated with the model for the mean signal count rate μs of  (2.12) is negligible. Now, we will relax these critical assumptions and approach the (perhaps) more realistic case when we must incorporate both the model error and the effects of an uncertain mean background count rate in the measurement model.Towards this objective, the following model is assumed for the (true though unknown) mean signal count rateμsj (j=1,2,…,m). The mean signal count rate can be described by the radiation data z and by the quantitative model (cf.  (2.12)), so
(2.13)μsj=μ̂sj+ej(1)=μ̃sj+ej(2),j=1,2,…,m,
where μ̂sj is an estimate for the mean signal count rate at the jth sensor obtained from the data z (more specifically, zj), μ̃sj is the mean signal count rate at the jth sensor obtained from the model given in (2.12), ej(1) is the uncertainty in the estimate μ̂sj of the mean signal count rate, and ej(2) represents the model error incurred by using (2.12) to predict the mean signal count rate. The model errors consist of a number of contributions which include (1) the effects of absorption of the gamma radiation by the air and (2) the effects of backscattering from the ground surface which may (possibly) influence the inverse-square law dependence in the model of (2.12).The measurement model of  (2.13) can be rearranged to give(2.14)μ̂sj-μ̃sj=ej(2)-ej(1)≡ϵj,j=1,2,…,m,
where ϵj is the composite error. If the expectation value 〈ϵj〉 of the composite error is zero and the variance 〈ϵj2〉 is assumed to be given by σj2, then application of the principle of maximum entropy (see, Jaynes [12]) provides the following explicit form for the likelihood function (assuming that at least one radiological source is present, so r>0):(2.15)p(z∣Θ,I)=1∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}.
It is noted that the Gaussian likelihood (which is effectively the probability distribution of the “noise” ϵj) obtained from application of the maximum entropy principle is simply the most conservative noise distribution that is consistent with the given facts (namely, the first two moments of the “noise” are known). In the Gaussian likelihood function of (2.15), the variance of the composite error ϵj is obtained as(2.16)σj2=σ1,j2+σ2,j2,j=1,2,…,m,
where σ1,j2 is the variance of the uncertainty ej(1) in the estimate μ̂sj of the mean signal count rate (obtained from the data z) and σ2,j2 is the variance of the model error ej(2). In this paper, the model error standard deviation (or square root of the variance) is specified as σ2,j=(μ̃sj)1/2.To complete the specification for the likelihood function in  (2.15), we need to provide the estimate μ̂sj along with the variance σ1,j2 in this estimate. To this purpose, we assume that we have available a measurement of background radiation (obtained in the absence of any point radiation sources) given to us as zb counts obtained over an exposure time τb. Now, suppose that a radiation sensor at a fixed receptor location obtains a measurement of z counts in an exposure time τ. We will now show how to use this information to obtain an estimate for the mean signal count rate μ̂s, as well as to characterize the uncertainty in this estimate in terms of the standard deviation σ1.To accomplish this, we utilize a Bayesian solution for the extraction of weak signals in strong background interference described by Loredo [13]. We will briefly outline the application of this methodology to the problem in this paper. To begin, we can apply Bayes' rule to compute the joint probability for the mean signal count rate μs and the mean background count rate μb conditioned on the measurement z. This gives(2.17)p(μs,μb∣z,I)∝p(μs∣μb,I)p(μb∣I)p(z∣μs,μb,I).
It is stressed that the contextual information I in (2.17) includes the contextual information Ib concerning the measurement of the background radiation and the contextual information Is concerning the measurement of the radiation field which may include the presence of one or more radiation sources. Hence, the contextual information I that is available in the problem can be decomposed as I=Ib,Is.We will assign a least informative prior probability forμs (the so-called Jeffreys prior; see Jaynes [12]) which has the following form (conditioned on a known mean background rate μb):(2.18)p(μs∣μb,I)=1μs+μb.
The prior probability for μb is chosen to be the following informative prior (using the information that is available, namely, that a measurement of the background radiation yielded zb counts in an exposure time τb):(2.19)p(μb∣I)≡p(μb∣zb,Ib,Is)=p(μb∣zb,Ib)∝p(μb∣Ib)p(zb∣μb,Ib),
since Is does not have any bearing on the determination of μb. The prior p(μb∣Ib) is assigned the nonsinformative Jeffreys prior so p(μb∣Ib)=1/μb and the likelihood function p(zb∣μb,Ib)=(μbτb)zbexp(-μbτb)/zb!, so(2.20)p(μb∣I)=p(μb∣zb,Ib)∝τb(μbτb)zb-1exp⁡⁡(-μbτb).The likelihood function in (2.17), associated with the measurement of z counts in exposure interval τ coming from a radiation field with mean count rate (μs+μb), is described by the following Poisson distribution:(2.21)p(z∣μs,μb,I)=τz(μs+μb)zz!exp⁡⁡(-(μs+μb)τ).
Now, if we substitute (2.18), (2.20), and (2.21) in  (2.17), the joint probability for the mean signal and background count rates is given by(2.22)p(μs,μb∣z,I)∝(μs+μb)z-1μbzb-1exp⁡⁡(-μsτ)exp⁡⁡(-μb(τb+τ)).
The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the (unknown) mean background count rate:(2.23)p(μs∣z,I)=∫0∞p(μs,μb∣z,I)dμb,
which on using  (2.22) and performing the integration finally gives (on normalization of the probability density function)(2.24)p(μs∣z,I)=∑i=1zΦiτ(μsτ)i-1(i-1)!exp⁡⁡(-μsτ),
where(2.25)Φi=(1+τb/τ)i(z+zb-i-1)!/(z-i)!∑j=1z(1+τb/τ)j(z+zb-j-1)!/(z-j)!.
Note that  (2.24) and (2.25) provide an explicit expression for the posterior probability of the mean signal count rate μs that is independent of the (unknown) mean background count rate μb. More specifically, the precise value of μb is not known, so it is treated as a nuisance parameter here and removed by integration (marginalization).Given the posterior probability ofμs summarized by (2.24) and (2.25), a best estimate for the mean signal count rate can be obtained as the posterior mean:(2.26)μ̂s≡〈μs〉=∫0∞μsp(μs∣z,I)dμs=1τ∑i=1ziΦi.
The second moment of μs about zero can be evaluated as(2.27)〈μs2〉=∫0∞μs2p(μs∣z,I)dμs=1τ2∑i=1zi(i+1)Φi,
from which the uncertainty σ1 in the estimate for μs given in  (2.26) can be determined as(2.28)σ1=(〈μs2〉-(μ̂s)2)1/2.
Equations (2.26), (2.27) and (2.28) when used in conjunction with (2.16) and (2.15) fully determine the form of the likelihood function (assuming implicitly that r>0). This functional form for the likelihood function accounts explicitly for the model error and the sampling error (the latter of which is compounded by an uncertain mean background count rate).For the caser=0, the likelihood function needs to be determined as follows (again assuming that our state of knowledge consists of a background radiation measurement yielding zb counts in exposure time τb):(2.29)p(z∣zb,I)=∏j=1m∫0∞P(zj∣μbτ)p(μb∣zb,Ib)dμb,
which on insertion of the normalized form of  (2.20) for p(μb∣zb,Ib) and evaluation of the integral gives(2.30)p(z∣zb,I)=∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!.
In (2.29) and (2.30), the radiation measurements zj are assumed to be made using the same exposure time τ (without any loss in generality).Finally, the likelihood function of (2.15) valid for r>0 and the likelihood function of (2.30) valid for r=0 can be combined to give the following general likelihood function for the measurement model assumed in this section:(2.31)p(z∣zb,Θ,I)=I{r>0}(r)∏j=1m2πσjexp⁡⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}+I{r=0}(r)∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!,
where the assumption of the availability of a background radiation count zb (obtained over an exposure time τb) has been included in the conditioning of the likelihood function (and of the posterior distribution based upon this likelihood function).


## 3. Computational Framework

This section describes the computational procedures that were used for extracting the source parameter estimates for source reconstruction. To this purpose, two different methodologies for sampling from two different forms of the posterior distribution for the source parameters are described.


### 3.1. Importance Sampling Using Progressive Correction

Importance sampling using progressive correction is a computational methodology for sampling from a posterior distribution. In this paper, this sampling procedure is applied to a form of the posterior distribution for the source parameters in which the number of sourcesr and the mean background count rate μb are assumed to be known a priori (see  (5.1)). To this purpose, the information contained in the prior distribution for the source parameters p(θ(r)∣r,I) and in the measurements z assimilated through the likelihood function p(z∣μb,r,θ(r),I)≡l(z∣μb,r,θ(r)) is combined to obtain the posterior distribution for θ(r) (assuming that the number of sources r and the background count rate μb are known):(3.1)p(θ(r)∣μb,r,z,I)∝l(z∣μb,r,θ(r))p(θ(r)∣r,I).
The minimum mean-squared error (MMSE) estimate of θ is then the posterior expectation,(3.2)θ(r)̂≡〈θ(r)∣μb,r,z〉=∫θ(r)p(θ(r)∣μb,r,z,I)dθ(r).Because the posterior PDFp(θ(r)∣μb,r,z,I) and, hence, the posterior expectation cannot be found exactly for the measurement model of Section 2, an approximation of the integral in (3.2) is computed via importance sampling. This involves drawing N samples of the parameter vector from an importance density and approximating the integral by a weighted sum of the samples:(3.3)θ(r)̂≈∑n=1NwnθS,n(r),
where {θS,n(r),n=1,…,N} is a “particle” set consisting of samples drawn from the importance density and {wn,n=1,…,N} is the set of importance weights. The sum of these weights is equal to one (namely, ∑n=1Nwn=1).Because the prior distribution will often be more diffuse than the likelihood, many samples drawn from the prior may fall outside the region of parameter space favoured by the likelihood function. Therefore, a straightforward application of (3.1) may result in poor source parameter estimates. The progressive correction approach overcomes this problem by using a multistage procedure in which samples are obtained from a series of posterior distributions which become progressively closer to the true posterior distribution. This is achieved by adopting an approximate likelihood at each stage which is somewhat more diffuse than the true likelihood.LetS denote the number of stages and ps(θ(r)∣μb,r,z,I), s=1,…,S denote the posterior PDF at stage s. A series of ps(θ(r)∣μb,r,z,I) can be constructed by setting(3.4)ps(θ(r)∣μb,r,z,I)∝ls(z∣μb,r,θ(r))p(θ(r)∣r,I),
where the intermediate likelihood at stage s=1,…,S is ls(z∣μb,r,θ(r))=l(z∣μb,r,θ(r))Gs with Gs=∑ℓ=1sγℓ, such that γs∈[0,1), for s=1,…,S and GS=α≤1. Assume that a random sample {wn,θs-1,n(r)} from ps-1(θ(r)∣μb,r,z,I) is available, and it is desired to produce a sample from ps(θ(r)∣μb,r,z,I). The steps of the PC algorithm are given in Algorithm 1. Note that for s=1, the sample {wn,θs-1,n(r)} is drawn directly from the prior p(θ(r)∣r,I) with wn=1/N.Algorithm 1:Importance sampling with progressive correction algorithm.
(1)s=0, G0=0, μb=μb*, r=r*;(2)forn=1,…,Ndo(3)   Drawθ0,n(r)~p(θ(r)∣r,I);(4)   Computel(z∣μb,r,θ0,n(r));(5)end for(6)whileGs<α and s<Sdo(7)s←s+1;(8)   Selectγs;(9)Gs=Gs-1+γs;(10)forn=1,…,Ndo(11)    Weights:ws,n=l(z∣μb,r,θs-1,n(r))γs/∑j=1Nl(z∣μb,r,θs-1,n(r))γs;(12)end for(13)  Resample{θs-1,n(r)}n=1N with weights {ws,n}n=1N to give {θn*}n=1N with uniform weights(14)forn=1,…,Ndo(15)θs,n(r)=θn*+ϵs,n, where ϵ~gs;(16)end for(17)end whileThe performance of the procedure depends somewhat on the number of stepsS and the expansion factors γ1,…,γS. The expansion factors γs were selected in line 8 of Algorithm 1 using an adaptive scheme proposed in Musso et al. [14]. In this adaptive scheme, the expansion factors are selected after each step rather than being selected a priori. Line 13 in Algorithm 1 performs a resampling on the samples {ws,n,θs-1,n(r)}n=1N to give samples {1/N,θn*}n=1N that have uniform weights. This resampling is undertaken using the systematic resampling algorithm (see, Kitagawa [15]), which is summarized in the form of a pseudocode given in Table  3.2 of Ristic et al. [16]. Line 15 in Algorithm 1 performs regularisation (jittering) of samples in order to avoid their duplication. Here, gs is a Gaussian regularisation kernel (see, Musso et al. [14]).The caseα=1 in the progressive correction algorithm should be used only when the likelihood function model in (2.9) is perfectly correct and its parameters precisely known. As we discussed in Section 2.2.2, this is unrealistic in practice: there are many potential errors in the model, such as the background radiation level, air attenuation, sensor location, and propagation loss. In order to make the PC algorithm robust against the modelling errors, it is recommended to adopt α<1 (this change affects the “while” loop in line 6 of Algorithm 1). In this way, the measurement likelihood is effectively approximated by a fuzzy membership function, but the Bayesian estimation framework holds as the weights in Step  9 of Algorithm 1 are normalised. This is similar to the treatment of the ambiguously generated unambiguous measurements in ([17], Chapter 7). Parameter α is a tuning parameter, which needs to be determined by field trials and calibration.The PC estimate of (3.2) is approximated as (after the resampling step of line 13 and the regularization step of line 15 in Algorithm 1)(3.5)θ(r)̂=1N∑n=1NθS,n(r).


#### 3.1.1. Source Number Estimation

The importance sampling algorithm using progressive correction described above assumes that the number of sources is knowna priori (see line 1 of Algorithm 1) which initializes r=r*, where r* is the number of sources. However, the number of sources will not be known in advance in practice, requiring this number also to be estimated from the data.For this purpose, we use the minimum description length (MDL) criterion which chooses the value ofr that minimises (Kay [18] page 225]); Rissanen [19]):(3.6)χr=-log⁡⁡l(z∣μb,r,θ̂ML(r))+nr2log⁡⁡m,
where θ̂ML(r) is the maximum likelihood estimate of θ(r) assuming r sources are present and nr≡3r is the dimension of θ(r) under the hypothesis r∈ℳ. Here, ℳ={0,…,rmax} denotes a set of candidate source numbers up to some maximum rmax. While MDL is derived based on the maximum likelihood estimate θ̂ML(r) (see  (3.6)), we will use the MDL in conjunction with the PC by simply replacing θ̂ML(r) with the PC estimate θ(r)̂ in  (3.6), so(3.7)χr=-log⁡⁡l(z∣μb,r,θ(r)̂)+nr2log⁡⁡m.
Finally, the estimate of the number of sources r is then determined as follows:(3.8)r̂=arg⁡min⁡r∈M⁡χr.To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values ofr=1,…,rmax. For each of these values of r, the estimated value for θ(r)̂ in  (3.5) is used to compute the corresponding value for χr from (3.7). The case r=0 (namely, the hypothesis that there are no sources present) is treated separately. There is no need to run the importance sampling using PC in this case, as (3.7) reduces to(3.9)χ0=-log⁡⁡l(z∣μb)
when r=0.


### 3.2. Reversible-Jump Markov Chain Monte Carlo Algorithm

In Section3.1, we described an importance sampling algorithm for sampling from the posterior distribution based on a sequence of proposal densities (involving, as such, a sequence of modified forms of a likelihood function that has been raised to powers γj<1, j=1,2,…,S). In its application here, calculations with the importance sampling algorithm using progressive correction were made separately for each model structure (namely, for each value of r), and comparisons of the different models were realized using the MDL criterion. In this section, we describe an alternative for sampling from the posterior distribution that involves application of an MCMC algorithm that allows sampling over both model and parameter spaces simultaneously (namely, r which indexes the various model structures is treated simply as another parameter, so that it is estimated jointly with the other model parameters θ(r)).Towards this objective, we apply a reversible-jump MCMC (RJMCMC) algorithm. The formalization of RJMCMC algorithms for dealing with variable dimension models has been described in the seminal work of Green [20]. The application of the RJMCMC algorithm to the problem of inverse dispersion has been described previously by Yee [21, 22]. The algorithm involves constructing a Markov chain whose stationary distribution is the posterior distribution p(Θ∣zb,z,I) of the source parameters Θ≡(rθ(r)⊺)⊺∈ℝ1+3r. To this purpose, let {Θ(t)}≡{(r(t)θ(r(t))⊺)⊺} (t=0,1,2,…) be the state vector of a Markov chain that is so constructed that its stationary distribution coincides with p(Θ∣zb,z,I) (see  (5.2) for the explicit form of the posterior distribution that is sampled from using the RJMCMC algorithm in this paper). The construction of the Markov chain uses an RJMCMC algorithm in which the MCMC moves over the model and parameter spaces are separated into two categories: (1) propagation moves which do not change the dimensionality of the parameter space; and, (2) transdimensional jump moves which change the source distribution by ±1 radiological point source and, in so doing, changes the parameter space dimension by ±3.In the first category of moves (which are dimension-conserving), the parameter vectorθ(r) (r fixed) is partitioned as follows: θ(r)=(θ1⊺θ2⊺)⊺ where θ1⊺≡(Qs,1Qs,2⋯Qs,r)∈ℝr and θ2⊺≡(xs,1⊺xs,2⊺⋯xs,r⊺)∈ℝ2r. The parameters collected together in θ1 are linearly related to the radiation “measurements” (cf.  (2.12)). For these parameters, we apply a Gibbs sampler for the update (propagation) move involving drawing samples directly from the univariate full conditional distribution for Qs,i (which happens to be a truncated Bernoulli-Gaussian distribution in this case). The remaining parameters in θ2 are related nonlinearly to the radiological measurements (cf.  (2.12)). Owing to the fact that the univariate conditional posterior distribution for these parameters cannot be determined analytically, we apply a Metropolis-Hastings (M-H) sampler to update the parameters in θ2. Towards this purpose, we update θk2∈θ2 (k=1,2,…,2r) by generating a new value θk2′ from a proposal transition kernel that is chosen here to be a Gaussian mixture, with each component of the mixture having a mean θk2 and different variances βk2. The variances of the components in this Gaussian mixture are chosen to cover several orders of magnitude from 0.001 to 0.1 times the “length” of the domain of definition for the parameter (assigned using the prior distribution for the parameter).In the second category of moves (which are dimension-changing), the source distribution is modified by±1 radiological point source by (1) a birth move that adds a single radiological point source and (2) a death move that removes a single radiological point source to/from the current source distribution. If a birth move is selected, the “coordinates” of the new radiological point source are generated by drawing a random sample from the prior distribution for each coordinate. Alternatively, if a death (reverse) move is selected, then a radiological point source in the current source distribution is randomly selected and removed. Following the recommendation of Green [20], the probabilities for the birth (pb) and death (pd) moves are specified as follows:(3.10)pb(r)=12min⁡⁡{1,p(r+1∣I)p(r∣I)},r=rmin⁡,…,rmax⁡-1,pd(r+1)=12min⁡⁡{1,p(r∣I)p(r+1∣I)},r=rmin⁡,…,rmax⁡-1,
ensuring that the probability of a jump move (either birth or death) pb+pd lies in [0.5,1] at each iteration. We note that for r=rmin, pd(r)=0 and for r=rmax, pb(r)=0.In practical terms, the RJMCMC algorithm as applied to the problem of source reconstruction can be set up as follows.(1)
Specify values for the hyperparameters(rmin,rmax,Qmax,γ,p*) and 𝒜 which define the prior distribution p(Θ∣I), as well as a value for tupper (maximum number of iteration steps to take).(2)
Set the iteration countert=1 and choose an initial state Θ(t-1) for the Markov chain by sampling from p(Θ∣I).(3)
Starting fromΘ(t-1), conduct the following sequence of moves to update the state vector to Θ(t):
(3.11)Θ(t-1)→Mb,dΘ⋆→M1Θ⋆⋆→M2Θ(t),
where Θ⋆ and Θ⋆⋆ denote some intermediate transition states between iterations (t-1) and t.(4)
Change the counter fromt to t+1, and return to step 3 until a maximum number of steps (tupper) has been taken.In Step  3,ℳb,d denotes a jump move involving either the birth (addition) of a radiological point source at a random locations or the death (removal) of an existing radiological point source with probabilities pb and pd, respectively; ℳ1 denotes the update of the source intensity rates using a Gibbs sampler; and, ℳ2 denotes the update of the source locations using an M-H sampler.


#### 3.2.1. Simulated Annealing

To improve convergence of the Markov chain, we have implemented a simulated annealing procedure. This procedure is very similar to the progressive correction that is used in conjunction with the importance sampling (as described in Section3.1). In a nutshell, simulated annealing mimics the quasistationary evolution of a system through a sequence of equilibrium samplings, involving an infinitely slow switching between the initial and final states. More concretely, let us consider an ensemble of N* different source distributions (system) that has been randomly drawn from a modified posterior having the following form: pλ(Θ∣zb,z,I)∝p(Θ∣I)pλ(z∣zb,Θ,I), where λ∈[0,1] is a “coolness” parameter. These samples will be labelled Θk(λ) (k=1,2,…,N*). When λ=0 (initial state), the likelihood function is switched off and the modified posterior distribution reduces exactly to the prior distribution. On the other hand, when λ=1 (final state) the modified posterior distribution is exactly the posterior that we wish to sample from. In between these two extremes, with λ∈(0,1), the effects of the radiation measurements z are introduced gradually through the modified (or “softened”) likelihood function pλ(z∣zb,Θ,I).The parameterλ can be interpreted as an inverse temperature parameter T (λ≡1/T), so λ∈[0,1] implies that T∈[1,∞]. The modified posterior pλ(Θ∣zb,z,I) (λ∈[0,1)) corresponds to “heating” the posterior distribution to a temperature T=1/λ>1. This heated distribution is flattened relative to the posterior distribution p(Θ∣zb,z,I), with the result that moves required to adequately cover the parameter space become more likely. In other words, this facilitates Markov chain mobility in the remote regions of the posterior distribution, ensuring a more reliable relaxation of the chain into the potentially (exponentially) small regions of the hypothesis space where the posterior distribution is large.We initialize the stochastic sampling scheme atλ=0 (infinite temperature) by randomly drawing N* samples of source distributions Θk(0) (k=1,2,…,N*) from the prior distribution p0(Θ∣zb,z,I). Given an ensemble of N* samples Θk(λ) that has achieved equilibrium (at temperature T=1/λ) with respect to the modified posterior pλ(Θ∣zb,z,I), an ensemble of N* samples Θk(λ+δλ) that is consistent with pλ+δλ(Θ∣zb,z,I) (at the reduced temperature T=1/(λ+δλ), δλ>0) can be obtained by using the weighted resampling method (see, Gamerman and Lopes [23]) applied to Θk(λ) (k=1,2,…,N*). To this purpose, each sample Θk(λ) is assigned a normalized importance weight wk as follows:
(3.12)wk=pλ+δλ(Θk(λ)∣zb,z,I)pλ(Θk(λ)∣zb,z,I)×(∑j=1N*wj)-1=pδλ(z∣zb,Θk(λ),I)×(∑j=1N*wj)-1,
for k=1,2,…,N*. We note that this procedure of computing weights is analogous to line 11 in Algorithm 1 for the importance sampling algorithm using progressive correction.Next, a resample is drawn from the discrete distribution concentrated at the samples{Θk(λ)}k=1N* with weights {wk}k=1N*. This resampling involves mapping the random discrete measure defined by {wk,Θk(λ)}k=1N* into the new random discrete measure defined by {1/N*,Θk*}k*=1N* having uniform weights (namely, Pr{Θk*=Θk(λ)}=wk with the result that Θk*~pλ+δλ(Θ∣zb,z,I) where ~ denotes “is distributed as”). The resampling conducted here is analogous to line 13 in Algorithm 1 for the importance sampling algorithm using progressive correction. The resampling used to produce the weighted replica ensemble acts as a selection mechanism in the sense that the number of replicas (copies) of a particular sample Θk(λ) is increased for those samples that are assigned a large weight wk (which, in turn, enhances the number of relevant configurations in the ensemble). Note that the weight determined according to  (3.12) provides a measure as to whether the state Θk(λ) contains new information about areas of the “energy” landscape that are favourable (reflected in larger values of the likelihood function).There are many ways to implement this resampling (with replacement) from{Θk(λ)}k=1N*, but in this paper we use the systematic resampling scheme described by Kitagawa [15] (namely, the same resampling scheme that is used in the PC algorithm described in Section 3.1). This scheme requires O(N*) time to execute and minimizes the Monte Carlo variation in the resample. For completeness, we briefly describe this resampling scheme as applied to the current problem. A vector (n1,n2,…,nN*) of copies of the samples Θk(λ) (k=1,2,…,N*) is obtained by computing a vector (ρ1,ρ2,…,ρN*) of the cumulative sums of N*×(w1,w2,…,wN*), generating a random draw from the uniform distribution u~𝒰([0,1]), and determining nk (k=1,2,…,N*) from(3.13)nk=⌊ρk+u⌋-⌊ρk-1+u⌋,k=2,3,…,N*-1,n1=⌊ρ1+u⌋,nN*=N*-⌊ρN*-1+u⌋,
where ⌊·⌋ denotes the “integer part of”. After resampling, the weights for each member of the resample are reset to wk=1/N* (namely, an equal weight is assigned to each member of the resample). After each resampling step to give an ensemble of N* samples Θk(λ+δλ) that is in equilibrium at temperature T=1/(λ+δλ) (approximately or better) with respect to pλ+δλ(Θ∣zb,z,I), we apply Nr (typically 10) Markov chain transitions Tλ+δλ to each member of the ensemble obtained from the resampling operation. These Markov chain transitions utilize the update scheme given by  (3.11) and leave pλ+δλ(Θ∣zb,z,I) invariant.It is necessary to specify the sequence ofδλ's used to move the ensemble of samples {Θk(λ)}k=1N* from λ=0 to λ=1. This sequence defines the annealing schedule for λ∈[0,1]. Taking a hint from the physical process of annealing, it is useful to allow the system to cool slowly. This allows the system to transition through a sequence of quasiequilibrium states and, in so doing, minimizes the probability of the system being trapped into a metastable state along the path of annealing. To mimic this process for the case of simulated annealing, at each step of the annealing schedule (starting from λ=0 where the stochastic sampling scheme is exploring the prior) the step size δλ is chosen so that(3.14)1N*∑k=1N*|wk-1|≈R,
where wk are the weights computed in accordance to  (3.12) and R is a small constant (R≪1). The constant R is chosen to provide a limiter on the maximum permissible size of δλ in each step of the annealing schedule. Note that unlike the PC algorithm, the sequence of λi∈[0,1] (i=0,1,2,…,K) with λi<λj for i<j and λ0=0, λK=1 that define the annealing schedule need not sum to α (namely, ∑i=0Kλi≠α in general).Whenλ=1, the annealing phase is complete (corresponding, as such, to the burn-in phase of the algorithm) and probabilistic exploration of the hypothesis space proceeds using the RJMCMC algorithm described at the end of Section 3.2. In the application of this algorithm following simulated annealing, Step  2 is modified as follows: instead of sampling Θ(0) from the prior distribution p(Θ∣I), these initial states are simply chosen to coincide with the samples from the ensemble of source distributions {Θk(λ=1)}k=1N* obtained at the end of the annealing phase.


## 3.1. Importance Sampling Using Progressive Correction

Importance sampling using progressive correction is a computational methodology for sampling from a posterior distribution. In this paper, this sampling procedure is applied to a form of the posterior distribution for the source parameters in which the number of sourcesr and the mean background count rate μb are assumed to be known a priori (see  (5.1)). To this purpose, the information contained in the prior distribution for the source parameters p(θ(r)∣r,I) and in the measurements z assimilated through the likelihood function p(z∣μb,r,θ(r),I)≡l(z∣μb,r,θ(r)) is combined to obtain the posterior distribution for θ(r) (assuming that the number of sources r and the background count rate μb are known):(3.1)p(θ(r)∣μb,r,z,I)∝l(z∣μb,r,θ(r))p(θ(r)∣r,I).
The minimum mean-squared error (MMSE) estimate of θ is then the posterior expectation,(3.2)θ(r)̂≡〈θ(r)∣μb,r,z〉=∫θ(r)p(θ(r)∣μb,r,z,I)dθ(r).Because the posterior PDFp(θ(r)∣μb,r,z,I) and, hence, the posterior expectation cannot be found exactly for the measurement model of Section 2, an approximation of the integral in (3.2) is computed via importance sampling. This involves drawing N samples of the parameter vector from an importance density and approximating the integral by a weighted sum of the samples:(3.3)θ(r)̂≈∑n=1NwnθS,n(r),
where {θS,n(r),n=1,…,N} is a “particle” set consisting of samples drawn from the importance density and {wn,n=1,…,N} is the set of importance weights. The sum of these weights is equal to one (namely, ∑n=1Nwn=1).Because the prior distribution will often be more diffuse than the likelihood, many samples drawn from the prior may fall outside the region of parameter space favoured by the likelihood function. Therefore, a straightforward application of (3.1) may result in poor source parameter estimates. The progressive correction approach overcomes this problem by using a multistage procedure in which samples are obtained from a series of posterior distributions which become progressively closer to the true posterior distribution. This is achieved by adopting an approximate likelihood at each stage which is somewhat more diffuse than the true likelihood.LetS denote the number of stages and ps(θ(r)∣μb,r,z,I), s=1,…,S denote the posterior PDF at stage s. A series of ps(θ(r)∣μb,r,z,I) can be constructed by setting(3.4)ps(θ(r)∣μb,r,z,I)∝ls(z∣μb,r,θ(r))p(θ(r)∣r,I),
where the intermediate likelihood at stage s=1,…,S is ls(z∣μb,r,θ(r))=l(z∣μb,r,θ(r))Gs with Gs=∑ℓ=1sγℓ, such that γs∈[0,1), for s=1,…,S and GS=α≤1. Assume that a random sample {wn,θs-1,n(r)} from ps-1(θ(r)∣μb,r,z,I) is available, and it is desired to produce a sample from ps(θ(r)∣μb,r,z,I). The steps of the PC algorithm are given in Algorithm 1. Note that for s=1, the sample {wn,θs-1,n(r)} is drawn directly from the prior p(θ(r)∣r,I) with wn=1/N.Algorithm 1:Importance sampling with progressive correction algorithm.
(1)s=0, G0=0, μb=μb*, r=r*;(2)forn=1,…,Ndo(3)   Drawθ0,n(r)~p(θ(r)∣r,I);(4)   Computel(z∣μb,r,θ0,n(r));(5)end for(6)whileGs<α and s<Sdo(7)s←s+1;(8)   Selectγs;(9)Gs=Gs-1+γs;(10)forn=1,…,Ndo(11)    Weights:ws,n=l(z∣μb,r,θs-1,n(r))γs/∑j=1Nl(z∣μb,r,θs-1,n(r))γs;(12)end for(13)  Resample{θs-1,n(r)}n=1N with weights {ws,n}n=1N to give {θn*}n=1N with uniform weights(14)forn=1,…,Ndo(15)θs,n(r)=θn*+ϵs,n, where ϵ~gs;(16)end for(17)end whileThe performance of the procedure depends somewhat on the number of stepsS and the expansion factors γ1,…,γS. The expansion factors γs were selected in line 8 of Algorithm 1 using an adaptive scheme proposed in Musso et al. [14]. In this adaptive scheme, the expansion factors are selected after each step rather than being selected a priori. Line 13 in Algorithm 1 performs a resampling on the samples {ws,n,θs-1,n(r)}n=1N to give samples {1/N,θn*}n=1N that have uniform weights. This resampling is undertaken using the systematic resampling algorithm (see, Kitagawa [15]), which is summarized in the form of a pseudocode given in Table  3.2 of Ristic et al. [16]. Line 15 in Algorithm 1 performs regularisation (jittering) of samples in order to avoid their duplication. Here, gs is a Gaussian regularisation kernel (see, Musso et al. [14]).The caseα=1 in the progressive correction algorithm should be used only when the likelihood function model in (2.9) is perfectly correct and its parameters precisely known. As we discussed in Section 2.2.2, this is unrealistic in practice: there are many potential errors in the model, such as the background radiation level, air attenuation, sensor location, and propagation loss. In order to make the PC algorithm robust against the modelling errors, it is recommended to adopt α<1 (this change affects the “while” loop in line 6 of Algorithm 1). In this way, the measurement likelihood is effectively approximated by a fuzzy membership function, but the Bayesian estimation framework holds as the weights in Step  9 of Algorithm 1 are normalised. This is similar to the treatment of the ambiguously generated unambiguous measurements in ([17], Chapter 7). Parameter α is a tuning parameter, which needs to be determined by field trials and calibration.The PC estimate of (3.2) is approximated as (after the resampling step of line 13 and the regularization step of line 15 in Algorithm 1)(3.5)θ(r)̂=1N∑n=1NθS,n(r).


### 3.1.1. Source Number Estimation

The importance sampling algorithm using progressive correction described above assumes that the number of sources is knowna priori (see line 1 of Algorithm 1) which initializes r=r*, where r* is the number of sources. However, the number of sources will not be known in advance in practice, requiring this number also to be estimated from the data.For this purpose, we use the minimum description length (MDL) criterion which chooses the value ofr that minimises (Kay [18] page 225]); Rissanen [19]):(3.6)χr=-log⁡⁡l(z∣μb,r,θ̂ML(r))+nr2log⁡⁡m,
where θ̂ML(r) is the maximum likelihood estimate of θ(r) assuming r sources are present and nr≡3r is the dimension of θ(r) under the hypothesis r∈ℳ. Here, ℳ={0,…,rmax} denotes a set of candidate source numbers up to some maximum rmax. While MDL is derived based on the maximum likelihood estimate θ̂ML(r) (see  (3.6)), we will use the MDL in conjunction with the PC by simply replacing θ̂ML(r) with the PC estimate θ(r)̂ in  (3.6), so(3.7)χr=-log⁡⁡l(z∣μb,r,θ(r)̂)+nr2log⁡⁡m.
Finally, the estimate of the number of sources r is then determined as follows:(3.8)r̂=arg⁡min⁡r∈M⁡χr.To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values ofr=1,…,rmax. For each of these values of r, the estimated value for θ(r)̂ in  (3.5) is used to compute the corresponding value for χr from (3.7). The case r=0 (namely, the hypothesis that there are no sources present) is treated separately. There is no need to run the importance sampling using PC in this case, as (3.7) reduces to(3.9)χ0=-log⁡⁡l(z∣μb)
when r=0.


## 3.1.1. Source Number Estimation

The importance sampling algorithm using progressive correction described above assumes that the number of sources is knowna priori (see line 1 of Algorithm 1) which initializes r=r*, where r* is the number of sources. However, the number of sources will not be known in advance in practice, requiring this number also to be estimated from the data.For this purpose, we use the minimum description length (MDL) criterion which chooses the value ofr that minimises (Kay [18] page 225]); Rissanen [19]):(3.6)χr=-log⁡⁡l(z∣μb,r,θ̂ML(r))+nr2log⁡⁡m,
where θ̂ML(r) is the maximum likelihood estimate of θ(r) assuming r sources are present and nr≡3r is the dimension of θ(r) under the hypothesis r∈ℳ. Here, ℳ={0,…,rmax} denotes a set of candidate source numbers up to some maximum rmax. While MDL is derived based on the maximum likelihood estimate θ̂ML(r) (see  (3.6)), we will use the MDL in conjunction with the PC by simply replacing θ̂ML(r) with the PC estimate θ(r)̂ in  (3.6), so(3.7)χr=-log⁡⁡l(z∣μb,r,θ(r)̂)+nr2log⁡⁡m.
Finally, the estimate of the number of sources r is then determined as follows:(3.8)r̂=arg⁡min⁡r∈M⁡χr.To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values ofr=1,…,rmax. For each of these values of r, the estimated value for θ(r)̂ in  (3.5) is used to compute the corresponding value for χr from (3.7). The case r=0 (namely, the hypothesis that there are no sources present) is treated separately. There is no need to run the importance sampling using PC in this case, as (3.7) reduces to(3.9)χ0=-log⁡⁡l(z∣μb)
when r=0.


## 3.2. Reversible-Jump Markov Chain Monte Carlo Algorithm

In Section3.1, we described an importance sampling algorithm for sampling from the posterior distribution based on a sequence of proposal densities (involving, as such, a sequence of modified forms of a likelihood function that has been raised to powers γj<1, j=1,2,…,S). In its application here, calculations with the importance sampling algorithm using progressive correction were made separately for each model structure (namely, for each value of r), and comparisons of the different models were realized using the MDL criterion. In this section, we describe an alternative for sampling from the posterior distribution that involves application of an MCMC algorithm that allows sampling over both model and parameter spaces simultaneously (namely, r which indexes the various model structures is treated simply as another parameter, so that it is estimated jointly with the other model parameters θ(r)).Towards this objective, we apply a reversible-jump MCMC (RJMCMC) algorithm. The formalization of RJMCMC algorithms for dealing with variable dimension models has been described in the seminal work of Green [20]. The application of the RJMCMC algorithm to the problem of inverse dispersion has been described previously by Yee [21, 22]. The algorithm involves constructing a Markov chain whose stationary distribution is the posterior distribution p(Θ∣zb,z,I) of the source parameters Θ≡(rθ(r)⊺)⊺∈ℝ1+3r. To this purpose, let {Θ(t)}≡{(r(t)θ(r(t))⊺)⊺} (t=0,1,2,…) be the state vector of a Markov chain that is so constructed that its stationary distribution coincides with p(Θ∣zb,z,I) (see  (5.2) for the explicit form of the posterior distribution that is sampled from using the RJMCMC algorithm in this paper). The construction of the Markov chain uses an RJMCMC algorithm in which the MCMC moves over the model and parameter spaces are separated into two categories: (1) propagation moves which do not change the dimensionality of the parameter space; and, (2) transdimensional jump moves which change the source distribution by ±1 radiological point source and, in so doing, changes the parameter space dimension by ±3.In the first category of moves (which are dimension-conserving), the parameter vectorθ(r) (r fixed) is partitioned as follows: θ(r)=(θ1⊺θ2⊺)⊺ where θ1⊺≡(Qs,1Qs,2⋯Qs,r)∈ℝr and θ2⊺≡(xs,1⊺xs,2⊺⋯xs,r⊺)∈ℝ2r. The parameters collected together in θ1 are linearly related to the radiation “measurements” (cf.  (2.12)). For these parameters, we apply a Gibbs sampler for the update (propagation) move involving drawing samples directly from the univariate full conditional distribution for Qs,i (which happens to be a truncated Bernoulli-Gaussian distribution in this case). The remaining parameters in θ2 are related nonlinearly to the radiological measurements (cf.  (2.12)). Owing to the fact that the univariate conditional posterior distribution for these parameters cannot be determined analytically, we apply a Metropolis-Hastings (M-H) sampler to update the parameters in θ2. Towards this purpose, we update θk2∈θ2 (k=1,2,…,2r) by generating a new value θk2′ from a proposal transition kernel that is chosen here to be a Gaussian mixture, with each component of the mixture having a mean θk2 and different variances βk2. The variances of the components in this Gaussian mixture are chosen to cover several orders of magnitude from 0.001 to 0.1 times the “length” of the domain of definition for the parameter (assigned using the prior distribution for the parameter).In the second category of moves (which are dimension-changing), the source distribution is modified by±1 radiological point source by (1) a birth move that adds a single radiological point source and (2) a death move that removes a single radiological point source to/from the current source distribution. If a birth move is selected, the “coordinates” of the new radiological point source are generated by drawing a random sample from the prior distribution for each coordinate. Alternatively, if a death (reverse) move is selected, then a radiological point source in the current source distribution is randomly selected and removed. Following the recommendation of Green [20], the probabilities for the birth (pb) and death (pd) moves are specified as follows:(3.10)pb(r)=12min⁡⁡{1,p(r+1∣I)p(r∣I)},r=rmin⁡,…,rmax⁡-1,pd(r+1)=12min⁡⁡{1,p(r∣I)p(r+1∣I)},r=rmin⁡,…,rmax⁡-1,
ensuring that the probability of a jump move (either birth or death) pb+pd lies in [0.5,1] at each iteration. We note that for r=rmin, pd(r)=0 and for r=rmax, pb(r)=0.In practical terms, the RJMCMC algorithm as applied to the problem of source reconstruction can be set up as follows.(1)
Specify values for the hyperparameters(rmin,rmax,Qmax,γ,p*) and 𝒜 which define the prior distribution p(Θ∣I), as well as a value for tupper (maximum number of iteration steps to take).(2)
Set the iteration countert=1 and choose an initial state Θ(t-1) for the Markov chain by sampling from p(Θ∣I).(3)
Starting fromΘ(t-1), conduct the following sequence of moves to update the state vector to Θ(t):
(3.11)Θ(t-1)→Mb,dΘ⋆→M1Θ⋆⋆→M2Θ(t),
where Θ⋆ and Θ⋆⋆ denote some intermediate transition states between iterations (t-1) and t.(4)
Change the counter fromt to t+1, and return to step 3 until a maximum number of steps (tupper) has been taken.In Step  3,ℳb,d denotes a jump move involving either the birth (addition) of a radiological point source at a random locations or the death (removal) of an existing radiological point source with probabilities pb and pd, respectively; ℳ1 denotes the update of the source intensity rates using a Gibbs sampler; and, ℳ2 denotes the update of the source locations using an M-H sampler.


### 3.2.1. Simulated Annealing

To improve convergence of the Markov chain, we have implemented a simulated annealing procedure. This procedure is very similar to the progressive correction that is used in conjunction with the importance sampling (as described in Section3.1). In a nutshell, simulated annealing mimics the quasistationary evolution of a system through a sequence of equilibrium samplings, involving an infinitely slow switching between the initial and final states. More concretely, let us consider an ensemble of N* different source distributions (system) that has been randomly drawn from a modified posterior having the following form: pλ(Θ∣zb,z,I)∝p(Θ∣I)pλ(z∣zb,Θ,I), where λ∈[0,1] is a “coolness” parameter. These samples will be labelled Θk(λ) (k=1,2,…,N*). When λ=0 (initial state), the likelihood function is switched off and the modified posterior distribution reduces exactly to the prior distribution. On the other hand, when λ=1 (final state) the modified posterior distribution is exactly the posterior that we wish to sample from. In between these two extremes, with λ∈(0,1), the effects of the radiation measurements z are introduced gradually through the modified (or “softened”) likelihood function pλ(z∣zb,Θ,I).The parameterλ can be interpreted as an inverse temperature parameter T (λ≡1/T), so λ∈[0,1] implies that T∈[1,∞]. The modified posterior pλ(Θ∣zb,z,I) (λ∈[0,1)) corresponds to “heating” the posterior distribution to a temperature T=1/λ>1. This heated distribution is flattened relative to the posterior distribution p(Θ∣zb,z,I), with the result that moves required to adequately cover the parameter space become more likely. In other words, this facilitates Markov chain mobility in the remote regions of the posterior distribution, ensuring a more reliable relaxation of the chain into the potentially (exponentially) small regions of the hypothesis space where the posterior distribution is large.We initialize the stochastic sampling scheme atλ=0 (infinite temperature) by randomly drawing N* samples of source distributions Θk(0) (k=1,2,…,N*) from the prior distribution p0(Θ∣zb,z,I). Given an ensemble of N* samples Θk(λ) that has achieved equilibrium (at temperature T=1/λ) with respect to the modified posterior pλ(Θ∣zb,z,I), an ensemble of N* samples Θk(λ+δλ) that is consistent with pλ+δλ(Θ∣zb,z,I) (at the reduced temperature T=1/(λ+δλ), δλ>0) can be obtained by using the weighted resampling method (see, Gamerman and Lopes [23]) applied to Θk(λ) (k=1,2,…,N*). To this purpose, each sample Θk(λ) is assigned a normalized importance weight wk as follows:
(3.12)wk=pλ+δλ(Θk(λ)∣zb,z,I)pλ(Θk(λ)∣zb,z,I)×(∑j=1N*wj)-1=pδλ(z∣zb,Θk(λ),I)×(∑j=1N*wj)-1,
for k=1,2,…,N*. We note that this procedure of computing weights is analogous to line 11 in Algorithm 1 for the importance sampling algorithm using progressive correction.Next, a resample is drawn from the discrete distribution concentrated at the samples{Θk(λ)}k=1N* with weights {wk}k=1N*. This resampling involves mapping the random discrete measure defined by {wk,Θk(λ)}k=1N* into the new random discrete measure defined by {1/N*,Θk*}k*=1N* having uniform weights (namely, Pr{Θk*=Θk(λ)}=wk with the result that Θk*~pλ+δλ(Θ∣zb,z,I) where ~ denotes “is distributed as”). The resampling conducted here is analogous to line 13 in Algorithm 1 for the importance sampling algorithm using progressive correction. The resampling used to produce the weighted replica ensemble acts as a selection mechanism in the sense that the number of replicas (copies) of a particular sample Θk(λ) is increased for those samples that are assigned a large weight wk (which, in turn, enhances the number of relevant configurations in the ensemble). Note that the weight determined according to  (3.12) provides a measure as to whether the state Θk(λ) contains new information about areas of the “energy” landscape that are favourable (reflected in larger values of the likelihood function).There are many ways to implement this resampling (with replacement) from{Θk(λ)}k=1N*, but in this paper we use the systematic resampling scheme described by Kitagawa [15] (namely, the same resampling scheme that is used in the PC algorithm described in Section 3.1). This scheme requires O(N*) time to execute and minimizes the Monte Carlo variation in the resample. For completeness, we briefly describe this resampling scheme as applied to the current problem. A vector (n1,n2,…,nN*) of copies of the samples Θk(λ) (k=1,2,…,N*) is obtained by computing a vector (ρ1,ρ2,…,ρN*) of the cumulative sums of N*×(w1,w2,…,wN*), generating a random draw from the uniform distribution u~𝒰([0,1]), and determining nk (k=1,2,…,N*) from(3.13)nk=⌊ρk+u⌋-⌊ρk-1+u⌋,k=2,3,…,N*-1,n1=⌊ρ1+u⌋,nN*=N*-⌊ρN*-1+u⌋,
where ⌊·⌋ denotes the “integer part of”. After resampling, the weights for each member of the resample are reset to wk=1/N* (namely, an equal weight is assigned to each member of the resample). After each resampling step to give an ensemble of N* samples Θk(λ+δλ) that is in equilibrium at temperature T=1/(λ+δλ) (approximately or better) with respect to pλ+δλ(Θ∣zb,z,I), we apply Nr (typically 10) Markov chain transitions Tλ+δλ to each member of the ensemble obtained from the resampling operation. These Markov chain transitions utilize the update scheme given by  (3.11) and leave pλ+δλ(Θ∣zb,z,I) invariant.It is necessary to specify the sequence ofδλ's used to move the ensemble of samples {Θk(λ)}k=1N* from λ=0 to λ=1. This sequence defines the annealing schedule for λ∈[0,1]. Taking a hint from the physical process of annealing, it is useful to allow the system to cool slowly. This allows the system to transition through a sequence of quasiequilibrium states and, in so doing, minimizes the probability of the system being trapped into a metastable state along the path of annealing. To mimic this process for the case of simulated annealing, at each step of the annealing schedule (starting from λ=0 where the stochastic sampling scheme is exploring the prior) the step size δλ is chosen so that(3.14)1N*∑k=1N*|wk-1|≈R,
where wk are the weights computed in accordance to  (3.12) and R is a small constant (R≪1). The constant R is chosen to provide a limiter on the maximum permissible size of δλ in each step of the annealing schedule. Note that unlike the PC algorithm, the sequence of λi∈[0,1] (i=0,1,2,…,K) with λi<λj for i<j and λ0=0, λK=1 that define the annealing schedule need not sum to α (namely, ∑i=0Kλi≠α in general).Whenλ=1, the annealing phase is complete (corresponding, as such, to the burn-in phase of the algorithm) and probabilistic exploration of the hypothesis space proceeds using the RJMCMC algorithm described at the end of Section 3.2. In the application of this algorithm following simulated annealing, Step  2 is modified as follows: instead of sampling Θ(0) from the prior distribution p(Θ∣I), these initial states are simply chosen to coincide with the samples from the ensemble of source distributions {Θk(λ=1)}k=1N* obtained at the end of the annealing phase.


## 3.2.1. Simulated Annealing

To improve convergence of the Markov chain, we have implemented a simulated annealing procedure. This procedure is very similar to the progressive correction that is used in conjunction with the importance sampling (as described in Section3.1). In a nutshell, simulated annealing mimics the quasistationary evolution of a system through a sequence of equilibrium samplings, involving an infinitely slow switching between the initial and final states. More concretely, let us consider an ensemble of N* different source distributions (system) that has been randomly drawn from a modified posterior having the following form: pλ(Θ∣zb,z,I)∝p(Θ∣I)pλ(z∣zb,Θ,I), where λ∈[0,1] is a “coolness” parameter. These samples will be labelled Θk(λ) (k=1,2,…,N*). When λ=0 (initial state), the likelihood function is switched off and the modified posterior distribution reduces exactly to the prior distribution. On the other hand, when λ=1 (final state) the modified posterior distribution is exactly the posterior that we wish to sample from. In between these two extremes, with λ∈(0,1), the effects of the radiation measurements z are introduced gradually through the modified (or “softened”) likelihood function pλ(z∣zb,Θ,I).The parameterλ can be interpreted as an inverse temperature parameter T (λ≡1/T), so λ∈[0,1] implies that T∈[1,∞]. The modified posterior pλ(Θ∣zb,z,I) (λ∈[0,1)) corresponds to “heating” the posterior distribution to a temperature T=1/λ>1. This heated distribution is flattened relative to the posterior distribution p(Θ∣zb,z,I), with the result that moves required to adequately cover the parameter space become more likely. In other words, this facilitates Markov chain mobility in the remote regions of the posterior distribution, ensuring a more reliable relaxation of the chain into the potentially (exponentially) small regions of the hypothesis space where the posterior distribution is large.We initialize the stochastic sampling scheme atλ=0 (infinite temperature) by randomly drawing N* samples of source distributions Θk(0) (k=1,2,…,N*) from the prior distribution p0(Θ∣zb,z,I). Given an ensemble of N* samples Θk(λ) that has achieved equilibrium (at temperature T=1/λ) with respect to the modified posterior pλ(Θ∣zb,z,I), an ensemble of N* samples Θk(λ+δλ) that is consistent with pλ+δλ(Θ∣zb,z,I) (at the reduced temperature T=1/(λ+δλ), δλ>0) can be obtained by using the weighted resampling method (see, Gamerman and Lopes [23]) applied to Θk(λ) (k=1,2,…,N*). To this purpose, each sample Θk(λ) is assigned a normalized importance weight wk as follows:
(3.12)wk=pλ+δλ(Θk(λ)∣zb,z,I)pλ(Θk(λ)∣zb,z,I)×(∑j=1N*wj)-1=pδλ(z∣zb,Θk(λ),I)×(∑j=1N*wj)-1,
for k=1,2,…,N*. We note that this procedure of computing weights is analogous to line 11 in Algorithm 1 for the importance sampling algorithm using progressive correction.Next, a resample is drawn from the discrete distribution concentrated at the samples{Θk(λ)}k=1N* with weights {wk}k=1N*. This resampling involves mapping the random discrete measure defined by {wk,Θk(λ)}k=1N* into the new random discrete measure defined by {1/N*,Θk*}k*=1N* having uniform weights (namely, Pr{Θk*=Θk(λ)}=wk with the result that Θk*~pλ+δλ(Θ∣zb,z,I) where ~ denotes “is distributed as”). The resampling conducted here is analogous to line 13 in Algorithm 1 for the importance sampling algorithm using progressive correction. The resampling used to produce the weighted replica ensemble acts as a selection mechanism in the sense that the number of replicas (copies) of a particular sample Θk(λ) is increased for those samples that are assigned a large weight wk (which, in turn, enhances the number of relevant configurations in the ensemble). Note that the weight determined according to  (3.12) provides a measure as to whether the state Θk(λ) contains new information about areas of the “energy” landscape that are favourable (reflected in larger values of the likelihood function).There are many ways to implement this resampling (with replacement) from{Θk(λ)}k=1N*, but in this paper we use the systematic resampling scheme described by Kitagawa [15] (namely, the same resampling scheme that is used in the PC algorithm described in Section 3.1). This scheme requires O(N*) time to execute and minimizes the Monte Carlo variation in the resample. For completeness, we briefly describe this resampling scheme as applied to the current problem. A vector (n1,n2,…,nN*) of copies of the samples Θk(λ) (k=1,2,…,N*) is obtained by computing a vector (ρ1,ρ2,…,ρN*) of the cumulative sums of N*×(w1,w2,…,wN*), generating a random draw from the uniform distribution u~𝒰([0,1]), and determining nk (k=1,2,…,N*) from(3.13)nk=⌊ρk+u⌋-⌊ρk-1+u⌋,k=2,3,…,N*-1,n1=⌊ρ1+u⌋,nN*=N*-⌊ρN*-1+u⌋,
where ⌊·⌋ denotes the “integer part of”. After resampling, the weights for each member of the resample are reset to wk=1/N* (namely, an equal weight is assigned to each member of the resample). After each resampling step to give an ensemble of N* samples Θk(λ+δλ) that is in equilibrium at temperature T=1/(λ+δλ) (approximately or better) with respect to pλ+δλ(Θ∣zb,z,I), we apply Nr (typically 10) Markov chain transitions Tλ+δλ to each member of the ensemble obtained from the resampling operation. These Markov chain transitions utilize the update scheme given by  (3.11) and leave pλ+δλ(Θ∣zb,z,I) invariant.It is necessary to specify the sequence ofδλ's used to move the ensemble of samples {Θk(λ)}k=1N* from λ=0 to λ=1. This sequence defines the annealing schedule for λ∈[0,1]. Taking a hint from the physical process of annealing, it is useful to allow the system to cool slowly. This allows the system to transition through a sequence of quasiequilibrium states and, in so doing, minimizes the probability of the system being trapped into a metastable state along the path of annealing. To mimic this process for the case of simulated annealing, at each step of the annealing schedule (starting from λ=0 where the stochastic sampling scheme is exploring the prior) the step size δλ is chosen so that(3.14)1N*∑k=1N*|wk-1|≈R,
where wk are the weights computed in accordance to  (3.12) and R is a small constant (R≪1). The constant R is chosen to provide a limiter on the maximum permissible size of δλ in each step of the annealing schedule. Note that unlike the PC algorithm, the sequence of λi∈[0,1] (i=0,1,2,…,K) with λi<λj for i<j and λ0=0, λK=1 that define the annealing schedule need not sum to α (namely, ∑i=0Kλi≠α in general).Whenλ=1, the annealing phase is complete (corresponding, as such, to the burn-in phase of the algorithm) and probabilistic exploration of the hypothesis space proceeds using the RJMCMC algorithm described at the end of Section 3.2. In the application of this algorithm following simulated annealing, Step  2 is modified as follows: instead of sampling Θ(0) from the prior distribution p(Θ∣I), these initial states are simply chosen to coincide with the samples from the ensemble of source distributions {Θk(λ=1)}k=1N* obtained at the end of the annealing phase.


## 4. Experimental Trial

A radiological field trial was conducted on a large, flat, and open area without obstacles within a military site in Puckapunyal, Victoria, Australia. An AN/PDR-77 radiation survey meter equipped with a gamma probe containing two Geiger-Müller tubes to cover both low and high ranges of dose rates was used to collect the radiation data. When connected to the AN/PDR-77 radiation survey meter, this gamma probe is capable of measuring gamma radiation dose rates from background to 9.99 Svh-1 without saturating [24], and it has a fairly flat response (~±10%) from 0.1 to above 1 MeV (see [24]; page 7). Measured dose rate data were recorded in μSv h-1. We converted these dose rate data into raw count measurements zj∈ℕ∪{0} by multiplication by a conversion factor (which was 0.21-1 for our probe) to convert to count rate data and subsequently by multiplication by the exposure time to convert to raw count data (measured over the exposure time).Three radiation sources were used in the field trial: two cesium sources (137Cs) and one cobalt (60Co) source (see Table 1). While the strengths of the radiation sources are typically characterised in terms of their activity in GBq (cf. Table 1), for simplicity we will characterise the strength (intensity rate) of a radiation point source by the expected dose rate at a distance of one meter from the source. The strengths of our three radiation sources calibrated in this manner were Qs=1912, 392.28, and 98.07 μSv m2h-1, for Sources 1, 2, and 3, respectively. We will use these values as the “true” source strengths in our source estimations.Table 1
Radiation sources used in the field trial.
SourceTypeActivity (MBq)1137Cs26×1032137Cs5×103360Co0.2×103The sources were mounted at the same height above the ground as the gamma probe. To ensure that radiation sources appear as isotropic in the horizontal plane, they were placed in a vertical configuration so that the handling rods were pointing up. Radiation dose measurements were collected on a grid of surveyed points that were carefully measured and marked beforehand in the local Cartesian coordinate system on the asphalt surface of the airfield. The data were acquired when the trolley-mounted gamma probe was positioned over individual grid points. During data collection at each grid point, the gamma probe remained stationary until sixty measurements were acquired. The exposure time for each radiation dose measurement (effectively the sampling interval) was kept constant at aboutΔ=0.8 s. Multiple measurements were collected at each location so that the candidate algorithms could be tested using several different batches of experimental data, or different total exposure intervals can be used in the source reconstruction by choosing the appropriate number of data points at each grid point (e.g., to test the source reconstruction algorithms for raw counts obtained over an exposure interval τ, we can simply select n* data points at each grid position such that τ≈n*Δ).Data sets were collected with one, two, and three radiation sources emplaced, respectively. These data sets are referred to as test sets 1, 2 and 3. The sources used when collecting these data sets and their locations in the local Cartesian coordinate system are listed in Table2.Table 2
The locations of radiation sources in the field trial.
Test setSource 1 LocationSource 2 LocationSource 3 Location1(11,10) m——2(11,10) m(3,50) m—3(11,10) m(3,50) m(41,5) mAn areal picture of the Puckapunyal airfield site where the field trial was conducted is shown in Figure1. The three red stars show the locations where the three radiation sources were emplaced. The white cross symbols indicate the grid points where data were collected. The x- and y-axis markings on the plot show the distances along these directions in meters (as well as the local Cartesian coordinate system that is used here to refer to locations of sources and grid positions).Figure 1
Aerial image of the Puckapunyal airfield site where the field trial was conducted. The red stars located at coordinates(11,10) m, (3,50) m, and (41,5) m of the local Cartesian coordinate system indicate where the three sources were emplaced for Test 3. The numbered grid points marked with crosses are the measurement points.In order to estimate the background radiation level in the field trial, measurements were collected at a number of grid positions in theabsence of radiation sources. At each of these grid positions, the detector was held stationary until 120 measurements were acquired. The background radiation counts were found to be spatially homogeneous and verify a Poisson distribution with mean background count rate of μb≈0.9 counts s-1.


## 5. Implementational Aspects

### 5.1. Importance Sampling Algorithm Using PC

The importance sampling algorithm using progressive correction was applied to the experimental field trial data for source reconstruction. Towards this purpose, this algorithm was used to sample from the following posterior distribution of the source parameters:(5.1)p(θ∣μb,r,z,I)∝p(θ∣r,I)p(z∣μb,r,θ,I)∝∏j=1mP(zj;λj(Θ))×∏i=1rQs,iκ-1exp⁡(-Qs,i/ψ)Γ(κ)ψκ×IA(xs,i)meas⁡(A),
which combines the prior distribution for source strengths given by (2.5) and the prior distribution for source location given by (2.7) with the likelihood function given by (2.9). The quantities λj(Θ) are determined in accordance withs (2.10). In this form of the posterior distribution, it is assumed that the mean background count rate μb and the number of sources r are known a priori.The importance sampling algorithm using progressive correction is applied to the experimental field trial data using the following values for the hyperparameters that define the prior distribution in (5.1). The gamma distribution that defines the prior for the source strengths has a shape parameter κ=1.5 and a scale parameter ψ=8000 μSv m2 h-1. This choice of the hyperparameters for the gamma distribution provides a prior for the source strengths that is broad enough to cover all likely source strength values. The prior for the source location is a uniform distribution over a specified rectangular-shaped area 𝒜, typically containing all measurement points. For the current application, the area 𝒜 is selected as the following rectangular region: 𝒜≡{(-50,200)⊗(-20,150)}, where ⊗ denotes Cartesian product.The PC algorithm is initialized by drawing a random sample in3r-dimensional parameter space from p(θ∣r,I). The number of samples in the PC algorithm is set to N=2500r.The expansion factors and the actual number of stages of the PC algorithm are data dependent. We have tuned the PC algorithm so that on average the number of stages is directly proportional to the number of sources. Parameterα was set to 1.


### 5.2. Reversible-Jump MCMC Algorithm

The reversible-jump MCMC algorithm (applied in conjunction with simulated annealing) was used to infer an unknown number of radiological point sources from the experimental field trial data. This involves drawing samples of radiological point source distributions from the following posterior distribution for the source parametersΘ:(5.2)p(Θ∣zb,z,I)∝p(Θ∣I)p(z∣zb,Θ,I)∝I{r>0}(r)∏j=1m2πσjexp⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}×(rmax⁡-rmin⁡)!(r-rmin⁡)!(rmax⁡-r)!p*(r-rmin⁡)(1-p*)rmax⁡-r×∏i=1r(1-γ)δ(Qs,i)+γI(0,Qmax⁡)(Qs,i)Qmax⁡×IA(xs,i)meas⁡(A)+I{r=0}(r)∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!×(1-p*)rmax⁡,which combines the likelihood function given by (2.31) with the prior for the number of sources given by (2.4), the prior for the source strength given by (2.6), and the prior for the source location given by (2.7). For the source reconstruction, the hyperparameters that define the prior distribution in (5.2) are assigned the following values: rmin=0, rmax=4, p*=0.25, γ=0.25, Qmax=75,000 μSv m2 h-1, and 𝒜={(-50,200)⊗(-20,150)} which is used to define the prior bounds for the location (xs,ys) of any radiological point source. Note that the prior bounds for the source location are exactly those assigned for the application of importance sampling using PC.The stochastic sampling algorithm was executed withN*=50 members of an ensemble of distributions of radiological point sources (originally, randomly drawn from the prior distribution p(Θ∣I)). During the burn-in (or annealing) phase of the algorithm, the control strategy used to determine the annealing schedule (see  (3.14)) was applied with R=0.03. After termination of the annealing at λ=1, a further tupper=1000 iterations of the RJMCMC procedure were applied to each of the N*=50 members of the ensemble of distributions of radiological point sources obtained at the end of the annealing phase. This gave 50,000 samples of distributions of radiological point sources drawn from the posterior distribution p(Θ∣zb,z,I) given in (5.2) (corresponding, as such, to the probabilistic exploration phase of the algorithm).The samples of radiological point source distributions drawn from the posterior distribution of (5.2) can be used to make various inferences on the values of the source parameters. For example, the number of sources can be obtained from the maximum a posteriori (MAP) estimate as follows:(5.3)r̂=arg⁡max⁡r⁡p(r∣zb,z,I),
where the (marginal) posterior probability for the number of sources is estimated from the samples as follows:(5.4)p̂(r∣zb,z,I)=1N′#{t:r(t)=r}=1N′∑t=1N′I{r}(r(t)).
Here, N′ is the number of samples Θ(1),Θ(2),… drawn from a sampled realization of the Markov chain. Similarly, using the Markov chain samples drawn from the posterior distribution, the posterior means of the various source parameters θi(r̂) (various components of θ(r̂); see  (2.2)) can be estimated by(5.5)θi(r̂)̂=∑t=1N′θr(t)i(t)I{r̂}(r(t))∑t=1N′I{r̂}(r(t)).
Similar estimates can be constructed for the posterior standard deviation of θi(r̂) which can be used as a measure of uncertainty in the recovery of this quantity. Alternatively, a p% credible (or highest posterior distribution (HPD)) interval that contains the source parameter θi(r̂) with p% probability, with the upper and lower bounds specified such that the probability density within the interval is everywhere larger than outside it, can be used as a measure of the uncertainty in the determination of this quantity.In addition to these statistical quantities, the RJMCMC algorithm applied in conjunction with simulated annealing permits the normalization constant (or evidence)Z for the posterior distribution of (5.2) to be determined2. This can be accomplished using an approach known as thermodynamic integration (see, von der Linden [25]).3 In this approach, Z can be evaluated as follows:(5.6)log⁡⁡Z=∫01〈log⁡⁡[p(z∣zb,Θ,I)]〉λdλ,
where 〈·〉 denotes the expectation operation taken with respect to the modified posterior pλ(Θ∣zb,z,I). This quantity can be computed during the simulated annealing phase of the stochastic simulation procedure. More specifically, the integral in (5.6) is discretized using the following Riemann sum:(5.7)log⁡⁡Z≈∑i=1K〈log⁡⁡[p(z∣zb,Θ,I)]〉λiδλi,
where δλi≡(λi-λi-1) is determined according to (3.14); and, the annealing schedule so defined provides the set of quadrature points (λ0=0,λ1,…,λK-1,λK=1).The evidenceZ can be used to determine the relative entropy (or, information gain) corresponding to the assimilation of the radiation measurements z, which involves updating our state of knowledge about the source embodied in the prior distribution to that embodied in the posterior distribution. To this end, the relative entropy H can be determined as follows (see, Cover and Thomas [26]): (5.8)H≡∫log⁡⁡(p(Θ∣zb,z,I)p(Θ∣I))p(Θ∣zb,z,I)dΘ=〈log⁡⁡[p(z∣zb,Θ,I)]〉-log⁡⁡(Z).
The relative entropy (of the posterior with respect to the prior) is a measure of the information gain provided by the receipt of the radiation measurements z. Geometrically, this can be visualized as the logarithm of the volumetric factor by which the prior has been compressed to give the posterior in the hypothesis space (the greater this compression, the greater the information gain embodied in z).


## 5.1. Importance Sampling Algorithm Using PC

The importance sampling algorithm using progressive correction was applied to the experimental field trial data for source reconstruction. Towards this purpose, this algorithm was used to sample from the following posterior distribution of the source parameters:(5.1)p(θ∣μb,r,z,I)∝p(θ∣r,I)p(z∣μb,r,θ,I)∝∏j=1mP(zj;λj(Θ))×∏i=1rQs,iκ-1exp⁡(-Qs,i/ψ)Γ(κ)ψκ×IA(xs,i)meas⁡(A),
which combines the prior distribution for source strengths given by (2.5) and the prior distribution for source location given by (2.7) with the likelihood function given by (2.9). The quantities λj(Θ) are determined in accordance withs (2.10). In this form of the posterior distribution, it is assumed that the mean background count rate μb and the number of sources r are known a priori.The importance sampling algorithm using progressive correction is applied to the experimental field trial data using the following values for the hyperparameters that define the prior distribution in (5.1). The gamma distribution that defines the prior for the source strengths has a shape parameter κ=1.5 and a scale parameter ψ=8000 μSv m2 h-1. This choice of the hyperparameters for the gamma distribution provides a prior for the source strengths that is broad enough to cover all likely source strength values. The prior for the source location is a uniform distribution over a specified rectangular-shaped area 𝒜, typically containing all measurement points. For the current application, the area 𝒜 is selected as the following rectangular region: 𝒜≡{(-50,200)⊗(-20,150)}, where ⊗ denotes Cartesian product.The PC algorithm is initialized by drawing a random sample in3r-dimensional parameter space from p(θ∣r,I). The number of samples in the PC algorithm is set to N=2500r.The expansion factors and the actual number of stages of the PC algorithm are data dependent. We have tuned the PC algorithm so that on average the number of stages is directly proportional to the number of sources. Parameterα was set to 1.


## 5.2. Reversible-Jump MCMC Algorithm

The reversible-jump MCMC algorithm (applied in conjunction with simulated annealing) was used to infer an unknown number of radiological point sources from the experimental field trial data. This involves drawing samples of radiological point source distributions from the following posterior distribution for the source parametersΘ:(5.2)p(Θ∣zb,z,I)∝p(Θ∣I)p(z∣zb,Θ,I)∝I{r>0}(r)∏j=1m2πσjexp⁡{-12∑j=1m(μ̂sj-μ̃sjσj)2}×(rmax⁡-rmin⁡)!(r-rmin⁡)!(rmax⁡-r)!p*(r-rmin⁡)(1-p*)rmax⁡-r×∏i=1r(1-γ)δ(Qs,i)+γI(0,Qmax⁡)(Qs,i)Qmax⁡×IA(xs,i)meas⁡(A)+I{r=0}(r)∏j=1mτbzbτzj(τ+τb)zj+zb(zj+zb-1)!zj!(zb-1)!×(1-p*)rmax⁡,which combines the likelihood function given by (2.31) with the prior for the number of sources given by (2.4), the prior for the source strength given by (2.6), and the prior for the source location given by (2.7). For the source reconstruction, the hyperparameters that define the prior distribution in (5.2) are assigned the following values: rmin=0, rmax=4, p*=0.25, γ=0.25, Qmax=75,000 μSv m2 h-1, and 𝒜={(-50,200)⊗(-20,150)} which is used to define the prior bounds for the location (xs,ys) of any radiological point source. Note that the prior bounds for the source location are exactly those assigned for the application of importance sampling using PC.The stochastic sampling algorithm was executed withN*=50 members of an ensemble of distributions of radiological point sources (originally, randomly drawn from the prior distribution p(Θ∣I)). During the burn-in (or annealing) phase of the algorithm, the control strategy used to determine the annealing schedule (see  (3.14)) was applied with R=0.03. After termination of the annealing at λ=1, a further tupper=1000 iterations of the RJMCMC procedure were applied to each of the N*=50 members of the ensemble of distributions of radiological point sources obtained at the end of the annealing phase. This gave 50,000 samples of distributions of radiological point sources drawn from the posterior distribution p(Θ∣zb,z,I) given in (5.2) (corresponding, as such, to the probabilistic exploration phase of the algorithm).The samples of radiological point source distributions drawn from the posterior distribution of (5.2) can be used to make various inferences on the values of the source parameters. For example, the number of sources can be obtained from the maximum a posteriori (MAP) estimate as follows:(5.3)r̂=arg⁡max⁡r⁡p(r∣zb,z,I),
where the (marginal) posterior probability for the number of sources is estimated from the samples as follows:(5.4)p̂(r∣zb,z,I)=1N′#{t:r(t)=r}=1N′∑t=1N′I{r}(r(t)).
Here, N′ is the number of samples Θ(1),Θ(2),… drawn from a sampled realization of the Markov chain. Similarly, using the Markov chain samples drawn from the posterior distribution, the posterior means of the various source parameters θi(r̂) (various components of θ(r̂); see  (2.2)) can be estimated by(5.5)θi(r̂)̂=∑t=1N′θr(t)i(t)I{r̂}(r(t))∑t=1N′I{r̂}(r(t)).
Similar estimates can be constructed for the posterior standard deviation of θi(r̂) which can be used as a measure of uncertainty in the recovery of this quantity. Alternatively, a p% credible (or highest posterior distribution (HPD)) interval that contains the source parameter θi(r̂) with p% probability, with the upper and lower bounds specified such that the probability density within the interval is everywhere larger than outside it, can be used as a measure of the uncertainty in the determination of this quantity.In addition to these statistical quantities, the RJMCMC algorithm applied in conjunction with simulated annealing permits the normalization constant (or evidence)Z for the posterior distribution of (5.2) to be determined2. This can be accomplished using an approach known as thermodynamic integration (see, von der Linden [25]).3 In this approach, Z can be evaluated as follows:(5.6)log⁡⁡Z=∫01〈log⁡⁡[p(z∣zb,Θ,I)]〉λdλ,
where 〈·〉 denotes the expectation operation taken with respect to the modified posterior pλ(Θ∣zb,z,I). This quantity can be computed during the simulated annealing phase of the stochastic simulation procedure. More specifically, the integral in (5.6) is discretized using the following Riemann sum:(5.7)log⁡⁡Z≈∑i=1K〈log⁡⁡[p(z∣zb,Θ,I)]〉λiδλi,
where δλi≡(λi-λi-1) is determined according to (3.14); and, the annealing schedule so defined provides the set of quadrature points (λ0=0,λ1,…,λK-1,λK=1).The evidenceZ can be used to determine the relative entropy (or, information gain) corresponding to the assimilation of the radiation measurements z, which involves updating our state of knowledge about the source embodied in the prior distribution to that embodied in the posterior distribution. To this end, the relative entropy H can be determined as follows (see, Cover and Thomas [26]): (5.8)H≡∫log⁡⁡(p(Θ∣zb,z,I)p(Θ∣I))p(Θ∣zb,z,I)dΘ=〈log⁡⁡[p(z∣zb,Θ,I)]〉-log⁡⁡(Z).
The relative entropy (of the posterior with respect to the prior) is a measure of the information gain provided by the receipt of the radiation measurements z. Geometrically, this can be visualized as the logarithm of the volumetric factor by which the prior has been compressed to give the posterior in the hypothesis space (the greater this compression, the greater the information gain embodied in z).


## 6. Results

To test the usefulness of the two approaches for source reconstruction, we have applied them to test data sets 1, 2, and 3 corresponding to one-, two-, and three-source examples, respectively. However, owing to space limitations, we will show only the results of application of the two approaches for test case 3. This test is arguably the most difficult case, but the results from the test are representative also for the other two test cases.The results for test case 3 will be shown for two different applications: namely, the radiation measurements at the grid positions (see Figure1) were obtained over exposure intervals of τ=16 and 48 s. In the application of importance sampling using PC, the mean background count rate μb was assumed to be exactly known (μb=0.9 counts s-1 was used in the algorithm). For the application of RJMCMC used in conjunction with simulated annealing, the first 24 s (τb=24 s) of a measurement of the background counts zb at receptor location (xl,yl)=(-40,-15) m (or, at grid position 1 as shown in Figure 1) was used in the source recovery.Firstly, we show the results of the two applications of the importance sampling algorithm using PC for test case 3. For each of these two applications, the algorithm (summarized in Algorithm1) was executed with fixed values of r=r* with r*=1, 2, and 3. Essentially, the assumption in the algorithm is that the number of sources is exactly known a priori, with r=r*. At the completion of this process, the MDL was computed according to (3.7) for r∈{0,1,2,3} and a “best” estimate for the number of sources was found using (3.8). In both applications of the importance sampling algorithm using PC here (namely, for τ=16 and 48 s), the minimum value of χr was found to correspond to r̂=3, an estimate that coincides with the correct number of sources for test case 3.The results of the source recovery for the two applications of the importance sampling algorithm are summarized in Table3 for r=3 (assumed correct number of sources). Here, the posterior mean, posterior standard deviation, and lower and upper bounds for the 95% HPD interval of the source parameters (xs,ys,Qs) are shown for τ=16 s (upper portion of Table 3) and for τ=48 s (lower portion of Table 3). In addition, histograms of the source parameters (xs,ys,Qs) for the three sources are exhibited in Figures 2 and 3, respectively, for τ=16 and 48 s.Table 3
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parametersxs,i (m), ys,i (m), and Qs,i (μSv m2 h-1) for i=1, 2, and 3 calculated from samples of distributions of radiological point sources drawn from the posterior distribution of  θ(r)(r=3) using importance sampling with progressive correction. The upper and lower parts of the table show the results of the source recovery using radiological measurements with exposure times of τ=16 and 48 s, respectively.
ParameterMeanStandard deviation95% HPDActualτ=16 si=1xs10.70.2(10.3,11.0)11.0ys9.90.2(9.5,10.4)10.0Qs193043(1848,2014)1912i=2xs2.10.6(1.0,3.2)3.0ys49.80.3(49.3,50.3)50.0Qs33031(269,392)392i=3xs40.90.1(40.7,41.1)41.0ys5.00.3(4.5,5.6)5.0Qs92.915.2(62.3,122.7)98.1τ=48 si=1xs10.90.01(10.89,10.90)11.0ys9.60.01(9.62,9.63)10.0Qs19061.3(1905,1907)1912i=2xs2.30.01(2.36,2.38)3.0ys49.60.01(49.60,49.62)50.0Qs3360.6(335.3,336.2)392i=3xs40.80.00(40.87,40.87)41.0ys5.30.01(5.26,5.28)5.0Qs83.60.73(83.0,84.2)98.1Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distributionp(θ∣μb,r,z,I) for τ=16 s and r=3 using the importance sampling algorithm with PC correction. (a, b, c) Histograms for the three parameters, namely xs-coordinate of source, ys-coordinate of source, and source intensity rate Qs that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.
(a)
Source 1(b)
Source 2(c)
Source 3Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distributionp(θ∣μb,r,z,I) for τ=48 s and r=3 using the importance sampling algorithm with PC correction. (a, b, c) Histograms for the three parameters, namely, xs-coordinate of source, ys-coordinate of source, and source intensity rate Qs that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.
(a)
Source 1(b)
Source 2(c)
Source 3For the application withτ=16 s, the algorithm successfully recovered the true parameters for all three sources to within the stated errors (e.g., the true parameters were all contained within a three-standard deviation interval about the best estimates of the parameters, and/or the true parameters were contained within the 95% HPD interval). However, this is not the case for the application of the algorithm with τ=48 s. More specifically, in this case it is seen that neither the three-standard deviation interval about the best estimates (based on the posterior mean) nor the 95% HPD interval for (xs,ys,Qs) for all three sources enclose the true values for these parameters. The reason for this is that the likelihood function in (5.1) that is used in the importance sampling algorithm with PC does not incorporate model error (namely, the model for the radiation measurements given by (2.12) is assumed to be exact). In consequence, for τ=48 s, the sampling error is smaller than the model error, so the latter becomes the dominant component in the error specification. The neglect of the model error in this case leads to underestimates for the precision at which some of the source parameters are recovered by the importance sampling algorithm.Next, we show the results for source reconstruction for test case 3 using the RJMCMC algorithm (applied in conjunction with simulated annealing) for the two applications withτ=16 and 48 s. Figure 4(a) displays a trace plot for the number of radiological point sources r in source distribution models sampled from the posterior distribution of Θ (for τ=48 s). The comparable result for τ=16 s is very similar (and hence, not shown). The trace plot here corresponds to the first 500 iterations of the RJMCMC sampler obtained during the probabilistic exploration phase (post-burn-in iterations following the completion of the annealing phase). Note that death moves for source distribution models, involving transitions from r=3 to 2 (or from r=2 to 1 and from r=1 to 0) do not occur. Furthermore, birth moves involving transitions from r=2 to 3 do not occur, either. On the other hand, dimension-changing moves (births or deaths) involving transitions from r=3 to 4 (and vice versa) are seen to occur. The residence time for occupation of a state with r=3 is seen to be significantly longer than that for a state with r=4. This is confirmed in Figure 4(b) which exhibits the probability distribution p(r)≡p(r∣zb,z,I) for the number of sources (estimated as a histogram from the 50,000 post-burn-in samples in accordance with  (5.4)).Trace plot (a) of the number of radiological point sourcesr in the first 500 samples drawn from p(Θ∣zb,z,I) during the probabilistic exploration phase using RJMCMC (after the completion of the annealing phase) and the corresponding posterior distribution (b) for the number of sources, p(r)≡p(r∣zb,z,I).
(a)(b)Figure4(b) shows that the most probable number of sources (MAP estimate) for this test case is 3 (r̂=3), which corresponds to the correct number of sources for this test. More specifically, r̂=3 is favored with a probability of about 0.7. Interestingly, it is noted that many of the distributions having four radiological point sources (which from Figure 4(b) are favored with probability of about 0.3) have one of the point sources “turned” off (namely, the source strength Qs for this source is identically zero). In fact, these distributions of radiological point sources should be considered as being identical to r=3 sources (and not r=4 sources owing to the fact that Qs=0 for one of the sources implies that that source is not present, at least as far as the detection of that source by a radiological sensor is concerned). If this classification for the number of radiological sources is used instead, then r̂=3 is now favored with a probability greater than about 0.995.Given the fact that our best estimate for the number of sources isr̂=3, we can now estimate the source parameters (xs,ys,Qs) corresponding to each of these three radiological point sources. Towards this purpose, we extracted all samples of distributions having exactly three sources (including all four-source samples, with one of the sources “turned off”). It is noted that an important issue in the interpretation of the results is the identifiability problem, which arises owing to the fact that the posterior distribution of Θ (cf.  (5.2)) is invariant under a relabelling of the identifiers used for each of the radiological point sources in the distribution. In consequence, we relabel the radiological point sources so that sources 1, 2, and 3 correspond to the (arbitrary) labelling of sources used in Table 2. After this re-labelling, Table 4 summarizes the posterior mean and standard deviation, as well as the lower and upper bounds for the 95% HPD interval of the source parameters, for each of the three radiological point sources for radiological measurements with exposure times τ=16 and 48 s.Table 4
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parametersxs,i (m), ys,i (m), and Qs,i (μSv m2 h-1) for i=1, 2, and 3 calculated from samples of distributions of radiological point sources drawn from the posterior distribution of Θ using the RJMCMC algorithm. The upper and lower parts of the table show the results of the source recovery using radiological measurements with an exposure times of τ=16 and 48 s, respectively.
ParameterMeanStandard deviation95% HPDActualτ=16 si=1xs8.70.8(7.2,10.6)11.0ys10.60.8(9.1,12.1)10.0Qs1404152(1109,1704)1912i=2xs3.00.6(3.0,5.0)3.0ys49.80.6(49.0,50.9)50.0Qs38581(228,548)392i=3xs40.90.5(39.9,41.6)41.0ys5.00.9(3.5,6.4)5.0Qs19446(105,286)98.1τ=48 si=1xs9.30.9(7.6,11.0)11.0ys10.20.8(8.7,11.9)10.0Qs1433156(1130,1742)1912i=2xs3.80.6(3.0,4.9)3.0ys49.80.6(49.0,50.9)50.0Qs38280(225,543)392i=3xs41.10.5(40.0,41.8)41.0ys4.80.8(3.5,6.3)5.0Qs20452(104,309)98.1From Table4, it is seen that the estimated values for the locations of the three sources agree well with the true locations of the sources in the reconstruction undertaken for both τ=16 and 48 s. Comparing these estimated values for the source locations with the actual source locations, we see also that the algorithm has adequately recovered the source locations to within the errors as represented by the 95% HPD interval for τ=16 and 48 s. Furthermore, the intensity rates of the sources are estimated quite well for both τ=16 and 48 s, although the accuracy in the recovery of these parameters is not as good as for the source location. Although the 95% HPD interval does not contain the true estimate for Qs for sources 1 and 3 for both τ=16 and 48 s, it should be noted that the 99% HPD interval does enclose the true value of Qs for these sources.An examination of Table4 shows that the recovery of the source parameters for τ=16 and 48 s using the RJMCMC algorithm gives results that are very similar, both in terms of their accuracy and their precision. This is in contrast with the results obtained using the importance sampling algorithm where it was found that there was an underestimate in the precision of recovery of the source parameters for τ=48 s. This difference in performance arises because the importance sampling algorithm was applied to the simpler posterior distribution of  (5.1) which does not incorporate the effects of model error, whereas the RJMCMC algorithm was applied to the posterior distribution of  (5.2) where the effects of model error are explicitly incorporated.A comparison of Tables3 and 4 shows that the source parameters are generally recovered with greater accuracy and better precision using the importance sampling algorithm than those recovered using the RJMCMC algorithm. The reason for this is because the importance sampling algorithm uses a posterior distribution (cf.  (5.1)) that is much more informative than the posterior distribution (cf.  (5.2)) used by the RJMCMC algorithm. In particular, the number of sources is assumed known a priori in the former, but is an unknown parameter in the latter; the mean background rate is assumed to be known in the former, but is unknown in the latter (where it is assumed that only a background count measurement zb taken over an exposure time τb is available), and the effects of model error on the source recovery are taken into account in the latter, but not in the former.Histograms of the source parameters (approximating the marginal probability distributions for the parameters), associated with each of the three identified radiological point sources, are exhibited in Figure5 for τ=48 s. Again, the comparable results for τ=16 s are not shown owing to the fact that they are very similar. A perusal of this figure suggests that the localisation of the three sources in both the xs- and ys-directions is generally very good. The recovery of the source intensity rates Qs is quite good, although both the accuracy and precision in this recovery are poorer than the recovery of the source locations.Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distributionp(Θ∣zb,z,I) for τ = 48 s using the RJMCMC algorithm. (a, b, c) Histograms for the three parameters, namely, xs-coordinate of source, ys-coordinate of source, and source intensity rate Qs that characterize sources 1, 2, and 3. In each frame, the solid vertical line indicates the true value of the parameter and the dashed vertical line corresponds to the best estimate of the parameter obtained as the posterior mean of the marginal posterior distribution for the parameter.
(a)
Source 1(b)
Source 2(c)
Source 3Finally, the information gainH (see (5.8)) obtained from the radiation measurements z for τ=16 and 48 s (for test case 3) was found to be H=53.4 and 60.1 natural units (nits), respectively.4 This implies that the information contained in the radiation measurements with exposure times τ=16 and 48 s allowed the “posterior volume” of the hypothesis space (volume of hypothesis space of reasonably large plausibility after receipt of the radiation data) to decrease by a factor of exp(H)≈1.55×1023 and 1.26×1026, respectively, relative to the “prior volume” of the hypothesis space (volume of the hypothesis space of reasonably large plausibility before the receipt of the radiation data). It is seen that using radiation data measured with exposure time of τ=48 s provided about 7 nits more information gain than that using radiation data measured with an exposure time of τ=16 s.


## 7. Conclusions

This paper presented two approaches for estimating an unknown number of point radiation sources and their parameters employing radiation measurements collected using a gamma radiation detector. The first approach used importance sampling with PC to sample the posterior distribution, and the second approach used the reversible-jump MCMC technique in conjunction with simulated annealing for this purpose. The two approaches were compared by applying them to experimental data collected in a recent field trial.It was demonstrated that both approaches perform well on the experimental data: in general, the number of sources is correctly determined and the parameters (e.g., location and intensity rates) that characterize each radiological point source are reliably estimated, along with a determination of the uncertainty (e.g., standard deviation, credible intervals) in the inference of the source parameters. However, it was shown that application of the first approach withα=1 did not adequately characterize the uncertainty in the source parameter estimates for the radiation measurements obtained with an exposure time τ=48 s. The defect here was not due to the progressive correction algorithm per se, but rather the result of the fact that the simple form of the posterior distribution of θ used with this algorithm assumed that the model error was negligible. This turned out not to be the case, and application of the RJMCMC with a more sophisticated form of the posterior distribution of Θ that incorporated explicitly the model error was shown to correct this defect.The results of this paper suggest that Bayesian probability theory is a powerful tool for the formulation of methodologies for source reconstruction. It should be stated that the application of importance sampling using PC to the posterior distribution of  (5.1) is fully Bayesian in character (assuming that r is known a priori), as is the application of the RJMCMC algorithm to the posterior distribution of  (5.2). However, the use of MDL in conjunction with the importance sampling algorithm to determine an estimate of the number of sources is not consistent with a strict Bayesian formalism.From this perspective, it should be noted that if the Bayes factor is used instead to determine the number of sources (namely, to determine which model of a distribution of radiological point sources to use), then the methodology would be fully consistent with a Bayesian formalism (see, Morelande and Ristic [8]). In fact, as discussed in Jaynes [12], application of the Bayes factor for model selection automatically incorporates Ockham's razor (namely, selecting the simplest model for a distribution of radiological point sources that “best” fits the available radiation data).Importance sampling using PC may also be carried out over the joint space ofr and θ in the same way as the RJMCMC to estimate the number of sources and parameters associated with these sources simultaneously and in a manner that is fully consistent with the Bayesian formalism (see Ristic et al. [27]).



---

*Source: 189735-2011-10-25.xml*