Comparison of Two Approaches for Detection and Estimation of Radioactive Sources

1 Physical Protection Section, Defence Research and Development Canada-Suffield, P.O. Box 4000 StnMain, Medicine Hat, AB, Canada T1A 8K6 2 Human Protection and Performance Division (HPPD), Defence Science and Technology Organization, 506 Lorimer Street, Fishermans Bend, VIC 3207, Australia 3 Intelligence Surveillance and Reconnaissance Division (ISRD), Defence Science and Technology Organization, 506 Lorimer Street, Fishermans Bend, VIC 3207, Australia


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 unknown a 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 Section 2. Section 3 describes the two approaches we use for radio-logical 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.

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: 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.

Assignment of Prior Probability
Here, we assume that a number r 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 location x s,i ≡ x s,i , y s,i in a Cartesian coordinate system with x s,i , y s,i ∈ A, where A ⊂ R 2 is some large region that is assumed to contain all the sources, ii its intensity rate source strength , Q s,i .
The parameter vector of source i is thus given by where denotes the matrix transpose.The source parameter vectors are collected into a stacked vector θ 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: Now, we need to assign explicit functional forms for each of the component prior probabilities in 2.3 .The prior distribution of r, p r | I , is assigned as a binomial distribution: r r min , r min 1, . . ., r max .In 2.4 , p * ∈ 0, 1 is the binomial rate, and, r min and r max are, respectively, the minimum and maximum number of sources.
Two alternative functional forms are used for the prior distribution of Q s,i i 1, 2, . . ., r .The first of the two alternatives assigns the prior for the source strength Q s,i to be a gamma distribution: 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 Q s,i is assigned a Bernoulli-uniform mixture model: where γ is the probability that the source is active Pr{Q s,i > 0} γ , δ x is the Dirac delta function, and Q max is an a priori upper bound on the expected source intensity rate.In 2.6 , I A x denotes the indicator characteristic function for set A, with Finally, the prior distribution for the source location x s,i is assigned to be uniformly distributed over the region A ⊂ R 2 that is assumed a priori to contain the source: where meas A is the area of the region A.

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 ∈ N ∪ {0} counts N ∪ {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 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.Let z j ∈ N ∪ {0} j 1, . . ., m be a measured count from the GM counter, taken at receptor location x l,j , y l,j ∈ R 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.

Measurement Model: Known Background Rate
The likelihood function is simply the joint density of the measurement vector z ≡ z 1 ••• z m 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: where λ j Θ is the mean radiation count for the jth sensor location: 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 μ j s at the jth sensor due to the radioactive point sources only is modelled in 2.10 as

2.12
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 μ j s μ j s .

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 μ j s 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 where μ j s is an estimate for the mean signal count rate at the jth sensor obtained from the data z more specifically, z j , μ j s is the mean signal count rate at the jth sensor obtained from the model given in 2.12 , e 1 j is the uncertainty in the estimate μ j s of the mean signal count rate, and e 2 j 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 where j is the composite error.If the expectation value j of the composite error is zero and the variance 2 j is assumed to be given by σ 2 j , 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 : 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 where σ 2 1,j is the variance of the uncertainty e 1 j in the estimate μ j s of the mean signal count rate obtained from the data z and σ 2 2,j is the variance of the model error e 2 j .In this paper, the model error standard deviation or square root of the variance is specified as σ 2,j μ j s To complete the specification for the likelihood function in 2.15 , we need to provide the estimate μ j s along with the variance σ 2 1,j 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 z b 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 It is stressed that the contextual information I in 2.17 includes the contextual information I b concerning the measurement of the background radiation and the contextual information I s 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 I b , I s .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 : 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 z b counts in an exposure time 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:

ISRN Applied Mathematics
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 The marginal probability distribution for the mean signal count rate can be derived by marginalizing with respect to the unknown mean background count rate: which on using 2.22 and performing the integration finally gives on normalization of the probability density function where 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: The second moment of μ s about zero can be evaluated as from which the uncertainty σ 1 in the estimate for μ s given in 2.26 can be determined as

2.28
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 case r 0, the likelihood function needs to be determined as follows again assuming that our state of knowledge consists of a background radiation measurement yielding z b counts in exposure time τ b : which on insertion of the normalized form of 2.20 for p μ b | z b , I b and evaluation of the integral gives In 2.29 and 2.30 , the radiation measurements z j 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: where the assumption of the availability of a background radiation count z b 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 .

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.

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 sources r 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 : The minimum mean-squared error MMSE estimate of θ is then the posterior expectation, Because the posterior PDF p θ 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: where {θ S,n r , n 1, . . ., N} is a "particle" set consisting of samples drawn from the importance density and {w n , n 1, . . ., N} is the set of importance weights.The sum of these weights is equal to one namely, N n 1 w n 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.
Let S denote the number of stages and p s θ r | μ b , r, z, I , s 1, . . ., S denote the posterior PDF at stage s.A series of p s θ r | μ b , r, z, I can be constructed by setting where the intermediate likelihood at stage s 1, . . ., S is r, z, I is available, and it is desired to produce a sample from p s θ r | μ b , r, z, I .The steps of the PC algorithm are given in Algorithm 1.Note that for s 1, the sample {w n , θ s−1,n r } is drawn directly from the prior p θ r | r, I with w n 1/N.
The performance of the procedure depends somewhat on the number of steps S and the expansion factors γ 1 , . . ., γ S .The expansion factors γ s were selected in line 8 of Algorithm 1 Compute l z | μ b , r, θ 0,n r ; 5 end for 6 while G s < α and s < S do 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 {w s,n , θ s−1,n r } N n 1 to give samples {1/N, θ n * } N n 1 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, g s 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 θ S,n r . 3.5

Source Number Estimation
The importance sampling algorithm using progressive correction described above assumes that the number of sources is known a priori see line 1 of Algorithm 1 which initializes ISRN Applied Mathematics 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 of r that minimises Kay 18 page 225 ; Rissanen 19 : where θ ML r is the maximum likelihood estimate of θ r assuming r sources are present and n r ≡ 3r is the dimension of θ r under the hypothesis r ∈ M. Here, M {0, . . ., r max } denotes a set of candidate source numbers up to some maximum r max .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 Finally, the estimate of the number of sources r is then determined as follows: To estimate the number of sources using the MDL criterion, we first run the importance sampling using progressive correction for all values of r 1, . . ., r max .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 when r 0.

Reversible-Jump Markov Chain Monte Carlo Algorithm
In Section 3.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 Θ | z b , z, I of the source parameters Θ ≡ r θ 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 Θ | z b , 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 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 Q s,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 θ 2  k ∈ θ 2 k 1, 2, . . ., 2r by generating a new value θ 2 k from a proposal transition kernel that is chosen here to be a Gaussian mixture, with each component of the mixture having a mean θ 2 k and different variances β 2 k .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 p b and death p d moves are specified as follows: ensuring that the probability of a jump move either birth or death p b p d lies in 0.5, 1 at each iteration.We note that for r r min , p d r 0 and for r r max , p b 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 r min , r max , Q max , γ, p * and A which define the prior distribution p Θ | I , as well as a value for t upper maximum number of iteration steps to take .
2 Set the iteration counter t 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 : where Θ and Θ denote some intermediate transition states between iterations t − 1 and t.
4 Change the counter from t to t 1, and return to step 3 until a maximum number of steps t upper has been taken.

In
Step 3, M 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 p b and p d , respectively; M 1 denotes the update of the source intensity rates using a Gibbs sampler; and, M 2 denotes the update of the source locations using an M-H sampler.

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 Section 3.1 .In a nutshell, simulat-ed 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: , 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 | z b , Θ, 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 λ Θ | z b , 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 Θ | z b , 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 p 0 Θ | z b , z, I .Given an ensemble of N * samples Θ k λ that has achieved equilibrium at temperature T 1/λ with respect to the modified posterior p λ Θ | z b , z, I , an ensemble of N * samples Θ k λ δλ that is consistent with p λ δλ Θ | z b , 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 w k as follows: 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 This resampling involves mapping the random discrete measure defined by {w k , Θ k λ } N * k 1 into the new random discrete measure defined by 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 w k 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 λ } N * k 1 , but in this paper we use the systematic resampling scheme described by Kitaga-wa 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 n 1 , n 2 , . . ., n N * 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 * × w 1 , w 2 , . . ., w N * , generating a random draw from the uniform distribution u ∼ U 0, 1 , and determining where • denotes the "integer part of".After resampling, the weights for each member of the resample are reset to w k 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 ISRN Applied Mathematics p λ δλ Θ | z b , z, I , we apply N r 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 λ δλ Θ | z b , z, I invariant.It is necessary to specify the sequence of δλ's used to move the ensemble of samples {Θ k λ } N * k 1 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 where w k 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, K i 0 λ 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 } N * k 1 obtained at the end of the annealing phase.

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 Sv h −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 z j ∈ N ∪ {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 137 Cs and one cobalt 60 Co 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   and 3, respectively.We will use these values as the "true" source strengths in our source estimations.The 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 Table 2.
An areal picture of the Puckapunyal airfield site where the field trial was conducted is shown in Figure 1.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 .
In order to estimate the background radiation level in the field trial, measurements were collected at a number of grid positions in the absence 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 .

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: 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 m 2 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 A, typically containing all measurement points.For the current application, the area A is selected as the following rectangular region: A ≡ { −50, 200 ⊗ −20, 150 }, where ⊗ denotes Cartesian product.
The PC algorithm is initialized by drawing a random sample in 3r-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.

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 Θ: 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: r min 0, r max 4, p * 0.25, γ 0.25, Q max 75, 000 μSv m 2 h −1 , and A { −50, 200 ⊗ −20, 150 } which is used to define the prior bounds for the location x s , y s 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 with N * 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 t upper 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 Θ | z b , 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: where the marginal posterior probability for the number of sources is estimated from the samples as follows: 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, ISRN Applied Mathematics the posterior means of the various source parameters θ i r various components of θ r ; see 2.2 can be estimated by

5.5
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 determined 2 .This can be accomplished using an approach known as thermodynamic integration see, von der Linden 25 . 3In this approach, Z can be evaluated as follows: where • denotes the expectation operation taken with respect to the modified posterior p λ Θ | z b , 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: 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 evidence Z 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
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 .

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 Figure 1 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 us-ed 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 z b at receptor location x l , y l −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 Algorithm 1 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 Table 3 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 x s , y s , Q s 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 x s , y s , Q s for the three sources are exhibited in Figures 2 and 3, respectively, for τ 16 and 48 s.
For 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 x s , y s , Q s 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  Figure 4 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 Q s 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 Q s 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 is r 3, we can now estimate the source parameters x s , y s , Q s 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 dis-tribution.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.
From Table 4, 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 Q s 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 Q s for these sources.
An examination of Table 4 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 Tables 3 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 z b 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 Figure 5 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 x s -and y s -directions is generally very good.The recovery of the source intensity rates Q s is quite good, although both the accuracy and precision in this recovery are poorer than the recovery of the source locations.
Finally, the information gain H 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. 4This 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 × 10 23 and 1.26 × 10 26 , 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.

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 ISRN Applied Mathematics 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 of r 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 .
τ b : p μ b | I ≡ p μ b | z b , I b , I s p μ b | z b , I b ∝ p μ b | I b p z b | μ b , I b , 2.19 since I s does not have any bearing on the determination of μ b .The prior p μ b | I b is assigned the nonsinformative Jeffreys prior so p μ b | I b 1/μ b and the likelihood function p z b | μ b , I b μ b τ b

Figure 1 :
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.

Figure 2 :
Figure 2:Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution p θ | μ 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 x s -coordinate of source, y s -coordinate of source, and source intensity rate Q s 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.

Figure 3 :
Figure 3:Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution p θ | μ 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, x s -coordinate of source, y s -coordinate of source, and source intensity rate Q s 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.

Figure 4 :
Figure 4: Trace plot a of the number of radiological point sources r in the first 500 samples drawn from p Θ | z b , 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 | z b , z, I .

Figure 5 :
Figure5: Inference of the parameters of the three radiological point sources obtained from samples drawn from the posterior distribution p Θ | z b , z, I for τ 48 s using the RJMCMC algorithm.a, b, c Histograms for the three parameters, namely, xs-coordinate of source, y s -coordinate of source, and source intensity rate Q s 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.

:
Importance sampling with progressive correction algorithm.

Table 1 :
Radiation sources used in the field trial.

Table 2 :
The locations of radiation sources in the field trial.

Table 3 :
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters x s,i m , y s,i m , and Q s,i μSv m 2 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 τ

Table 4 :
The posterior mean, posterior standard deviation, and lower and upper bounds of the 95% HPD interval of the parameters x s,i m , y s,i m , and Q s,i μSv m 2 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.