The node-wise Pseudo-marginal method: model selection with spatial dependence on latent graphs

Motivated by problems from neuroimaging in which existing approaches make use of “mass univariate” analysis which neglects spatial structure entirely, but the full joint modelling of all quantities of interest is computationally infeasible, a novel method for incorporating spatial dependence within a (potentially large) family of model-selection problems is presented. Spatial dependence is encoded via a Markov random field model for which a variant of the pseudo-marginal Markov chain Monte Carlo algorithm is developed and extended by a further augmentation of the underlying state space. This approach allows the exploitation of existing unbiased marginal likelihood estimators used in settings in which spatial independence is normally assumed thereby facilitating the incorporation of spatial dependence using non-spatial estimates with minimal additional development effort. The proposed algorithm can be realistically used for analysis of moderately sized data sets such as 2D slices of whole 3D dynamic PET brain images or other regions of interest. Principled approximations of the proposed method, together with simple extensions based on the augmented spaces, are investigated and shown to provide similar results to the full pseudo-marginal method. Such approximations and extensions allow the improved performance obtained by incorporating spatial dependence to be obtained at negligible additional cost. An application to measured PET image data shows notable improvements in revealing underlying spatial structure when compared to current methods that assume spatial independence.


Introduction
Incorporating spatial dependence in the process of analysing large image, and other spatial image-like, data sets can be a difficult problem largely due to the computational requirements.An important example of such data is provided by PET (Positron Emission Tomography) imaging of the brain, where a whole image typically requires analysis of up to 10 6 time series (Hammers et al., 2007).Current state of the art analysis, e.g.Fan et al. (2021); Zhou et al. (2016), of such data generally either assumes spatial independence of pixels or performs large-scale aggregation over space (Zhou et al., 2002) to overcome computational restrictions.In this paper, a novel computational method to incorporate spatial dependence in the process of model selection is presented.
The technique presented here extends the well-known pseudo-marginal Markov Chain Monte Carlo (MCMC) algorithm originally presented by Beaumont (2003), and characterised by Andrieu and Roberts (2009).More specifically, a Potts model (Potts, 1952) is first utilised as a prior distribution to encode spatial dependence over a local model associated with each node of a graph.Imposing a further assumption of conditional spatial independence for local model parameters given the local model then gives rise to significant simplifications in computations.
This Markov random field model is used on the discrete space of model orders to describe spatial dependence between neighbouring nodes (pixels).Finally, a standard component-wise MCMC updating scheme is used in conjunction with the above assumptions.This then allows for the use of node-level marginal likelihood estimators to be used within arXiv:2109.08573v2[stat.ME] 17 Apr 2022 The NWPM Method a pseudo-marginal MCMC method.The end result is a tractable, accessible computational method that uses model selection at the individual (node/pixel) level to perform spatially dependent model selection on the whole (graph/image) data set.This method provides a flexible but efficient algorithm that can be readily implemented -in fact, the proposed method allows for the use of existing unbiased non-spatial estimators to be used within a framework which incorporates spatial dependence.
In addition, the proposed method can be further extended through augmentation of the state space of the MCMC chain; in this specialised setting careful specification of the proposal distributions can lead to many further computational approaches and techniques.
The remainder of this paper is structured as follows.Section 2 presents relevant background material including the motivating problem and reviews the use of Potts models in related contexts.Section 3 reviews existing methods for sampling from distributions of this type.Section 4 presents a hierarchical model and develops methods for using it in the context of interest alongside a formal justification of the approach and a number of extensions.Finally, the performance of the developed method, and some variants thereof, is studied empirically in applied settings in Section 5.

Background 2.1 Preliminaries
A graph G = (V, E) is the pair of sets of nodes, or vertices, denoted V and edges denoted E. In particular, each element of the edge-set E is some pair of elements of the nodes-set u, v ∈ V , denoted u, v ∈ E. Two nodes u, v are said to be connected by the edge u, v if u, v ∈ E. In this case we may say that u and v are neighbours, or adjacent, and denote this by the relation u ∼ v. Here, we will look only at undirected graphs, so the relation ∼ will be symmetric and u, v is an unordered pair.Denote by ∂(v) = {u ∈ V : v ∼ u} the set of neighbours of v, by convention v is not a neighbour of itself.
Given a graph, G, a collection of random variables, X = (X v : v ∈ V ), indexed by the nodes-set, V , is called a random field on G. Let P be the law of X; and define, for any A ⊂ V , With a slight abuse of this notation, we will write X −v to denote X V \{v} for v ∈ V .X is called a Markov random field (Besag, 1974) on a discrete graph if and only if we have that: where each component X v takes value in some set X .Thus, X takes values in X V , the collection of maps from V to X .In this paper, X will always be finite.
Denote a parametric model for data Y ∈ Y, by S = (Y, {P θ : θ ∈ Θ}); where P θ is a probability distribution over Y, with θ taking values in the parameter space Θ ⊂ R d with d < ∞.We will use f (•|θ) to denote the density (w.r.t some dominating measure, usually the Lebesgue or counting measures) of the distribution P θ .Allow Y ∼ P θ and, where the density exists, Y ∼ f (•|θ), to mean that Y is distributed according to the distribution P θ .
A finite set of statistical models, called the model space, is denoted S = {S M : M ∈ M}.Here, each model is indexed by the model order M , which is an element of a finite index set M ⊂ N.For a model S M , of order M , let f M (•|θ) for parameter θ ∈ Θ M , be the associated density.We write Θ M to emphasis that the parametric space is often dependent on the model order.A straightforward example of such models are the compartmental models, see (Gunn et al., 2001) or Section 2.2 below.Bayesian model selection, averaging and comparison involves inference of the model order M ∈ M, where M is a random variable to which some prior distribution is attached (Robert, 2007, Chapter 7).An important component of the Bayesian model selection process is the posterior model probability, denoted π(M |Y ).Specifically, a central objective here is to characterise this distribution in the setting for image and image-like data, while accounting for spatial dependence.The definition of the posterior distribution (particularly for the setting we are interested in) will be detailed in the sequel.

Position Emission Tomography
An application that motivates the methodological development that follows is Positron Emission Tomography (PET; Phelps et al. (2006)).
Dynamic PET is an increasingly important neuroimaging technique used to study different functions of the brain in vivo.PET images have been used to investigate many biochemical processes in the brain, in both control and atypical subjects.This invasive imaging modality uses the detection of high-energy photons, released as a consequence of an injected tracer's positron emissions, as a signal to form a dynamical three-dimensional image.
Given that there will be greater radioactivity in regions where there are higher amounts of the tracer, the signal intensity in each voxel will be proportional to the concentration of the tracer.Thus, by taking a sequence of measurements over time, the constructed image represents the time course of the tissue concentration of the tracer in vivo.
As a result of the nature and number of observations recorded in the reconstructed PET image, the data is usually modelled using a system of linear differential equations such as compartmental models (Gunn et al., 2001).We restrict our attention here to the plasma input compartment models.Given data y = (y 1 , . . ., y k ) , at a voxel, an M −compartment model can be formally written as: C T (t j ; φ 1:M , ϑ for j = 1, . . ., k; where C T is the tissue concentration, t j is the measurement time of observation y j , j is the additive measurement error and the plasma input function C P is treated as known.The parameters φ 1:M = (φ 1 , . . ., φ M ) and ϑ 1:M = (ϑ 1 , . . ., ϑ M ) determine the dynamics of the model.Importantly, these parameters are indistinguishable, (Gunn et al., 2001, Theorem 2.2), and are often called the micro-parameters.In contrast, the macro-parameter termed the volume of distribution is uniquely identifiable.This macro-parameter is very often the quantity of principal interest in statistical analysis of PET data.
The additive error variable j is typically modelled using a normal distribution: j ∼ N (0, σ 2 ), where N (0, σ 2 ) is the normal distribution with mean zero and variance σ 2 .This class of statistical models have proven to be very effective in Bayesian analysis of PET data, as studied and presented by Peng et al. (2008) and Jiang et al. (2009).
One ubiquitous method for PET analysis is the non-negative non-linear least squares (NNLS; Cunningham and Jones (1993)).This approach is often used due to the non-negative nature of the rate constants.The NNLS method involves minimal modelling assumptions, imposing no prior assumptions on the number of components.Instead, the time-series at each voxel is interpreted as a noisy measurement of exponential decay.The method involves the formulation of a constrained linear optimisation problem, where: M will represent the maximum number of terms to be included, for example Cunningham and Jones (1993) use M = 100; Each component of the ϑ parameter is fixed and predetermined to lie within some range that is physiologically meaningful.Next, the optimal value of φ j for each ϑ j is estimated using a numerical method, subject to φ i ≥ 0, j = 1, . . ., m.
Recently, computer-intensive Monte Carlo methods have also been successfully applied to these models towards meaningful statistical inference of PET data.Specifically, Zhou et al. (2013) investigated an application specific MCMC approach to Bayesian model comparison for this class of models; and studied both vague non-informative and biologically informed priors.In this work, samples from a MCMC chain were used in the generalized harmonic mean estimator to approximate the model evidence in order to perform model comparison, selection and averaging.This work also demonstrated the possibility that Monte Carlo approaches could be used to investigate and meaningfully compare more complex models, containing up to M = 3 compartments, rather than the typical up to M = 2 compartmental models.Zhou et al. (2013) further showed that within a Bayesian framework, a t−distributed error structure may be far more plausible.This can be written as j ∼ T (0, τ, ν), where T (0, τ, ν) denotes the Student's t−distribution with mean zero, scale τ , and ν degrees of freedom.Expressions for prior and posterior densities, of both cases, can be found in Appendix A.
Importantly, with regard to model selection within this context, Zhou et al. (2013) also showed that a Bayesian approach exhibited some spatial structure despite assuming spatial independence.This was in contrast to using the Akaike Information Criteria (AIC; Akaike (1973)) for model selection with a non-linear least square (NLS) method (as well as the NNLS method) used to approximate the MLE which showed no obvious spatial structure, see also Section 5.3.2.
More recently, Zhou et al. (2016) used an adaptive variant of the sequential Monte Carlo (SMC) sampler algorithm of Del Moral et al. (2006) to estimate the model evidence (marginal likelihood).The class of algorithms presented in this work, use adaptive MCMC kernels together with adaptive annealing schemes to minimise tuning requirements to further the accessibility of Monte Carlo approaches to this problem.In addition, the model evidence can be directly approximated using the set of weights from the sampler.The computational method proposed in this present paper, builds upon this approach and thus inherits many of these advantages.Further, Zhou et al. (2016) presented empirical results suggesting that this SMC method produces similar performance to the MCMC method proposed by Zhou et al. (2013); As such, in this paper, we focus on comparing performance to the SMC sampler method -which we refer to as the (spatially) independent SMC method.
Similarly, Castellaro et al. (2017) propose a method based on a Variational Bayes (VB) approach for parameter estimation in PET data analysis; and also present an empirical comparison to the MCMC method of Zhou et al. (2013).Specifically, in this empirical comparison the M = 2−compartment model was used to represent measured PET data.The VB method was adapted to kinetic modelling of PET data, applied to the M = 2−compartment model and used to estimate the V D .Next, the MCMC method, restricted to inference from the M = 2−compartment model, was also used to analyse this data.Castellaro et al. concluded that there was a strong agreement between the outputs of these two methods.In particular the V D maps of both the VB and MCMC methods showed a very high correlation (Pearson R 2 = 0.99 (see, Figure S2, Castellaro et al. (2017))).
There continues to be further development in Bayesian approaches in analysis of PET data, albeit in almost all cases they are non-spatial approaches.The most recent example is Fan et al. (2021), who propose a simple and intuitive algorithm, based on Approximate Bayesian Computation(ABC), for analysis of PET data.As with the above Bayesian approaches, this method, termed PET-ABC, assumes voxels are spatially independent.PET-ABC works by simulating from the parameter space using a simple rejection scheme: Firstly, proposals are sampled from the prior distribution, typically this is the uniform density with a physiologically meaningful range.Each proposal is used as a trial value to estimate the signal C T .Next, the error between a summary statistics of the estimated signal and the de-noised data y is computed.In their study, Fan et al. (2021) suggest using spline smoothed estimates of C T as the summary statistic for both the estimated and observed signal.The proposed value of the parameter is accepted if the above error is below a predetermined threshold.The generated parameter sample can then be use for point estimation and to quantify uncertainty.Bayesian model selection follows the above method, with the additional step of proposing a model at each iteration.Fan et al. (2021) state that using PET-ABC produced estimates with lower variance when compared to the weighted NLS method.

Potts Model
The Potts model (Potts, 1952), a generalisation of the Ising model (Ising, 1925), was used originally to model interacting spins on a lattice.However, these models have been shown to be also very effective in analysis of image data; for a detailed discussion and review of such applications see Winkler (1995), Geman (1990) and references therein.For example, Geman and Geman (1984) is an early study demonstrating the effectiveness of the Ising model for restoration of images under a Bayesian framework.More recently, the Potts model has been used very successfully within the more broader but still growing sub-field of Bayesian image analysis Hurn et al. (2003).
In the present context, the Potts model is an ideal prior for model orders, due to its discrete state space, minimal parametrisation and ability to encode the general principle that nearby vertices are a priori likely to be best described by the same model.Following the Bayesian approach, in this study, the model order of the data will be Markov random field with a Potts prior distribution.
For a review of Potts distributions, particularly in the image analysis context, see Hurn et al. (2003).
Given a graph G = (V, E), finite state space X and coupling constant J > 0, the Potts model specifies a family of parametric probability distributions on X V ; characterised by the joint mass function: Here, recall v ∼ u denotes neighbouring pairs and δ xv,xu is the Kronecker delta, i.e. δ xv,xu is one if x v = x u , and zero otherwise.The intractable normalising constant (or partition function), written The parameter J dictates how likely neighbouring random variable X v and X u are to have the same value.Given that the normalisation depends upon this parameter, J is typically very difficult to infer; and will be treated as known in this setting as discussed at the end of this section.
Given some graph G = (V, E), and the associated random variable X = (X v ∈ X : v ∈ V ), we write to mean that the random variable X, taking values in state space X V , is distributed according to the Potts model with given coupling constant J.Where no ambiguity arises, the parameters will be omitted from notation in the interests of clarity.
The Potts model exhibits phase transition behaviour.This has significant implications when using single-site update scheme MCMC methods, which could result in slow mixing.That is, if the chain is in a configuration such as the case described above, changes to single variables will mean proposing to move to states of lower probability.Furthermore, phase transition behaviour happens close to or higher than critical values of the coupling constant, which we denote as J critical .In the Monte Carlo setting, there is a sharp transition from fast mixing below J critical and slow mixing above it.Onsager (1944) showed that for the Ising model on the two dimensional first order square lattice the exact value was More recently, see for instance Matveev and Shrock (1996) or Wu (1982), this has been extended: In applications, since the intractable normalising constant ζ(J) is dependent on J, it is difficult to infer in the model fitting process (Everitt, 2012).Within the Monte Carlo context, these types of distributions are referred to as doublyintractable.Møller et al. (2006) introduced an ingenious method to address such problems for the Ising model.For a more recent work, and for a concise but detailed review of developments since Møller et al. (2006), see Moores et al. (2020) and references therein.
Additionally, J is a parameter of a prior distribution, thus under the Bayesian analysis framework it would not be too unreasonable for it to be pre-specified.As such, in what follows, we will treat J as known.For example, in the numerical studies below we use a value which gives good performance in all our pilot studies, with a preference for choosing a smaller value of J to avoid imposing too much spatial structure via this prior distribution.

Existing Monte Carlo Methods
Having established the problem of interest and reviewed the relevant distribution for encoding spatial information, we now turn our attention to a brief review of existing computational methods for investigating the characteristics of such complex distributions.We begin by looking at the relatively simple Gibbs sampling method for the Potts model.However, ultimately we would like to sample from, and subsequently characterise, a posterior distribution that consists of the Potts model as a prior, together with an intractable likelihood (detailed in Section 4).To this end, a common and standard approach to sampling from densities expressed in terms of intractable integrals, known as the pseudo-marginal method, is also discussed.In this paper, we use the SMC normalising constant estimator, to approximate said intractable integrals.Thus, the section ends with a brief discussion on the SMC estimator and its properties.However, generic pseudo-marginal methods would be computationally infeasible in the context of high-dimensional discrete latent spaces of the sort considered here; this then leads to the proposal of a novel computational method in the subsequent section.

Gibbs Sampler for the Potts Model
There exist efficient samplers for the Potts model, most notably Glauber Dynamics Glauber (1963) corresponding to single-site Gibbs updates, and the elegant Swendsen-Wang algorithm (Swendsen and Wang, 1987).Typically, the Swendsen-Wang algorithm is used in simple Potts models, where there is an absence of a "strong external filed" i.e. the (graph) likelihood in this context.In this regime, the Swendsen-Wang algorithm has appealing convergence properties for all values of the coupling strength, J, see Hurn et al. (2003).However, this excellent performance deteriorates markedly in the presence of an external field (Higdon, 1998) and hence is less appealing in the presence of an informative likelihood.In settings in which the likelihood is relatively uninformative, there might be benefits in combining this type of cluster update with the local moves explored here.
We focus on the Gibbs sampling approach (Glauber, 1963;Geman and Geman, 1984), variants of which will extend naturally to the context of interest in this paper.The full conditionals required for a Gibbs sampler are readily computed for the Potts model: As is typical in Gibbs samplers, the use of full conditionals results in a simple node-wise update schedule.
However, computing full conditionals for posterior densities such as that of interest here is not possible because the likelihood cannot be evaluated even up to a normalising constant.So we turn to another Monte Carlo method to tackle this problem.

Pseudo-Marginal Monte Carlo Markov Chain
Consider, first, MCMC in the simpler case of a tractable target density; one which can be evaluated point-wise.
Denote the target probability density µ(x), where x ∈ X , and X ⊆ R d is the state space.MCMC algorithms construct a µ−invariant Markov chain denoted (X (i) ) n i=1 and like all time-homogeneous Markov chains are completely characterised by their initial distribution and their transition kernel.
The Metropolis-Hastings (MH) (Metropolis et al., 1953;Hastings, 1970) algorithm is a particularly widely-used MCMC algorithm that may be used to target any analytically tractable density µ.The algorithm can be described very briefly as follows: Given that the MH Markov chain is at some state x ∈ X , a new state x * ∈ X is proposed from a proposal density, denoted here q(x, •).The proposed state x * is accepted with probability min{1, R(x, x * )}, where the acceptance ratio is defined R(x, x * ) := µ(x * )q(x * , x) µ(x)q(x, x * ) ; or remains in the current state otherwise.The MH algorithm specifies the kernel to be of the form: where, is the probability of remaining in the same state and δ x denotes a probability measure concentrated at x.
The requirement of µ being known point-wise arises when computing R. Thus, a very natural approach for intractable target densities is to approximate µ instead.Beaumont (2003) and Andrieu and Roberts (2009) show that as long as the approximation of µ is unbiased, the Markov chain of this approximate version will have the same invariant distribution as when the "exact" target distribution had been used.This approximation results in a class of techniques called "pseudo-marginal" algorithms.
Many intractable target densities, including those of interest here, can be expressed in terms of marginals of tractable densities: integrals that cannot be evaluated analytically.An obvious example is when we wish to use marginal likelihoods as will be the case in our context; see (10) below.
More formally, let the target density be of the form where the integrals cannot be solved analytically; and θ ∈ Θ may be considered a latent (or nuisance) variable that is not necessarily of current interest.
Let µ(x) be an estimator of µ(x) and suppose that µ(x) is unbiased, for all x ∈ X .That is, if we let g(•|x) denote the density of µ(x), then we have that E[ µ(x)] = µ(x) for every x ∈ X , where the expectation is taken with respect to g(•|x).Note that g need not be known.
The estimator µ(x) is random, and its variance will play an important part in producing accurate results -we postpone further discussion until Section 4.2.The normalising constant estimator of the sequential Monte Carlo sampler, as detailed in Section 3.3, is a common and popular marginal likelihood estimator.
Grouped independence Metropolis Hastings (GIMH), a pseudo-marginal algorithm first presented by Beaumont (2003) and subsequently interpreted by Andrieu and Roberts (2009), can now be presented.Using the simpler reformulation (without auxiliary variables) given by Andrieu and Vihola (2015), the pseudo-code description of GIMH is presented in Algorithm 1.
Here, we have used subscripts (i) and * for µ to emphasise precisely where the marginal likelihood has been (re-)estimated.Henceforth these subscripts will be suppressed for notational clarity.
2. Sample: 4. With probability min{1, r} let: Note in particular, Algorithm 1 is analogous to the marginal MH algorithm, with the difference being that the acceptance ratio, R, is now approximated by R(x, x * ) := µ(x * )q(x * , x) µ(x)q(x, x * ) .
As Andrieu and Roberts (2009) point out, the pseudo-marginal Markov chain should be intuitively thought of as a Markov chain on the extended space of the random variable (X, µ(X)), rather than just X itself.They further show that a Markov kernel with this acceptance ratio produces a Markov chain with (and converging to) invariant distribution with marginal µ.A similar justification applies to the algorithm proposed below which is a natural extension of GIMH specialised to the setting of interest.We detail this in Section 4.2.2,Proposition 1.

Estimating the marginal likelihood using an SMC Sampler
Unbiased likelihood estimators are commonly constructed using particle filters, or SMC samplers, as in Andrieu et al. (2010).
SMC samplers (Del Moral et al., 2006) are a class of algorithms that generate collections of weighted samples approximating each of a sequence of target distributions {µ t } t≥1 over essentially any random variable defined on some spaces (X t ) t≥1 .Standard SMC methods achieve this through a combination of sequential IS and variance reduction resampling methods on spaces with increasing dimensions (see, for example, Doucet and Johansen (2011)).However, SMC samplers allow for sequence of distributions defined on some common space X .
SMC samplers accomplish this by creating a sequence of auxiliary distributions on spaces of increasing dimensions.SMC samplers can also be readily applied to simpler spaces through the use of annealing schemes.See Zhou et al. (2016) and Section 5 for details on specification of the annealing scheme within the present Bayesian analysis context.The sequence of auxiliary distributions are specified using Markov kernels called "backward kernels", which critically impact the variance of the estimators.See Del Moral et al. (2006) for more details on these so-called backward kernels.
Conventional sequential importance re-sampling algorithms can then be applied to the auxiliary distribution sequence.More precisely, assume that at iteration t − 1 the set of N weighted particles approximate µ t−1 , denote this {W denotes the normalised weight of particle X (i) 0:t for all i = 1, . . ., N .For iteration t, the path of each particle X (i) 0:t−1 is extended with a Markov kernel say, K t (x t−1 , x t ), yielding the set of particles {X (i) 0:n } N i=1 to which IS is applied.Subsequently, the weights are updated by a factor called incremental weights.Importantly, particularly when interested in approximating marginals such as the model evidence (marginal likelihood); if µ t (x t ) = γ t (x t )/Z t is only known up to a normalising constant, the unnormalised incremental weights .
can be used.Note: here, and as practised in the studies below, we assume that K t is constructed to be µ t -invariant, µ t µ t−1 and the associated time-reversal kernel is used for the backward kernel.
The NWPM Method In fact, this collection of algorithms provide unbiased (see Del Moral (2004, Proposition 7.4.1))estimates of Z t /Z t−1 via These normalising constant estimators satisfy a Central Limit Theorem (Del Moral, 2004, Proposition 9.4.1) and have asymptotic variance that are inversely proportional to the number of particles N (see also Chopin and Papaspiliopoulos (2020) for a recent review).As such, N is an important tuning parameter when this sampler is used within a pseudomarginal algorithm (Doucet et al., 2015).

Methodology
We now turn our attention to the development and presentation of the models to encode spatial dependence together with computational methods that can be used to perform inference under, and characterise, these models.
There is a considerable literature on spatial modelling in Bayesian data analysis (see, for example, Gelfand et al. ( 2010)).However, some specialised settings, such as that encountered when dealing with PET images, are yet to be explored to the same extent.In these cases, there may be significantly large amounts of data to be analysed and relatively sophisticated models may be required for meaningful investigation.Such approaches may need to incorporate application-specific structure and information.For example, Bezener et al. (2018) presents a method of variable selection which models spatial dependence using hierarchical spatial priors and parcellation for the analysis of MRI images.
Studies that attempt to model spatial relationships in PET images are scarce, and tend not to use a Bayesian approach.
One such example is Zhou et al. (2002) who proposes a method to allow for spatial dependence.That is, spatial constraint using nonlinear ridge regression is used to improve upon parametric images produced by conventional weighted NLS methods.They showed that doing so reduced the percentage mean square error by 60 − 80% in simulation studies.These promising results motivate the incorporation of spatial dependence when analysing PET images, particularly in Bayesian frameworks where such studies have been largely absent.
On the other hand, there exists many computationally efficient methods for complex settings, such as PET data analysis.These include, in the PET context, among many others: Gunn et al. (2002); Peng et al. (2008); Jiang et al. (2009); Zhou et al. (2013Zhou et al. ( , 2016)).Computational tractability in these methods is typically attained through the assumption of spatial independence.
A natural strategy is, therefore, to adapt existing non-spatial methods to a model that does incorporate spatial dependence.This strategy is particularly natural when modelling spatial dependence via a Potts model as it can be efficiently targeted by MCMC methods with single-site update schemes.In this section we construct a class of models that allows for such computational methods.
Having presented the model that allows for effective incorporation of spatial dependence, we then present a natural extension of the generic pseudo-marginal approach, discussed in Section 3 above, to this specialised setting in Section 4.2 along with the theoretical justification of the method, which follows from that of the standard pseudo-marginal approach and some additional considerations.Further extensions of the method are described in Section 4.3, and a straightforward approximation of the proposed algorithm is introduced in Section 4.4.

Spatial Bayesian Model Selection using a Potts Prior
Given spatial data (or spatio-temporal data, such as PET For example, for image data a simple (finite) lattice would suffice.We will call Y v the node data point at the node v ∈ V .
Building good statistical models for Y can be difficult; model selection would involve selecting among many complex and computationally intensive models from the model space.Analysis in this setting can be simplified by using existing simpler models such as those that assume spatial independence and model data at the node-level.
Such an approach has two significant strengths.Firstly, this approach can readily exploit many existing techniques which make spatial-independence assumptions, allowing them to incorporate spatial dependence.Secondly, using this approach allows for considerable computational simplifications as will be seen later.

A Generic Model for Spatial Model Selection
Formally, consider parametric models for the node data point Y v , for all v ∈ V Y , rather than the whole data set Y .Denote the parameter of these models by θ v , with parametric space Θ v for each node v ∈ V .
It is immediate that a Potts model can be used to encode spatial relationships between models at different locations in space; define Here, the coupling constant J is treated as known and constant.
The parametric distribution for Y v , now has parameters (θ v , X v ).Here, X v is a component of the Potts random variable X, and so has spatial dependence dictated by E Y .The model order at node v is given by X v .Given a family of parameterised prior densities for each possible model denoted generically p, over the parameter space; we may compactly summarise this hierarchical model as: Here, recall that f M denotes the density associated with statistical model S M with model order M .A similar convention is used for the (parameter) prior to denote that a different prior distribution must be used for each model order in the Bayesian context.Specifically, here we have that the model order is M = X v for each node v ∈ V .This generic model is represented as a plate diagram in Figure 1.
Here, ξ v is a parameter of the distribution p.In principle a hyperprior could be attached to this parameter with no particular difficulty; we work with this simple setting in the interest of parsimony.To this end, we fix ξ v to be a known hyperparameter common to all nodes, v, and so we write ξ instead of ξ v in what follows.
Additionally, given that we now treat the model order as a random variable, henceforth we use the convention that Note that, at the nodal (pixel) level, Bayesian model selection would often involve the computation of the marginal likelihood, f ( We are interested in the marginal likelihood of the whole image Y , which may involve multiple intractable high-dimensional integrals.We impose an assumption of conditional independence, discussed next, which allows us to make significant simplifications.For a simple example, see the toy model described in Section 5.1; This toy model will be used for simulation experiments to evaluate the proposed method. In realistic settings, it may be that the parameter θ = (θ v ) v∈V also exhibits spatial dependence.The above hierarchical model, makes the underlying assumption that it does not and this setting is the focus of the present paper; where we have found the incorporation of spatial structure at the level of model sufficient to improve upon the state of the art for problems of interest.Further generalisation would be possible and provides an interesting direction for further work.We encode this in Assumption 1 which will be a standing assumption throughout this paper and is discussed in Section 4.1.3.
Assumption 1 (Conditional Independence) The primary aim of the proposed methodology is to infer X.Although θ is inferred as a by-product of the proposed method, it can be thought of as a latent variable in this framework and may be a nuisance parameter in some settings which mitigates the impact of modelling the parameters in this way.
Finally, note that, as a direct consequence of the above model and assumption, the marginal likelihood of Y can be written as the product of the marginal likelihoods of Y v over all v ∈ V -see Section 4.1.4below.

Spatial Bayesian Model Selection
To exploit the presented model, denote by S v = {S v,Mv : M v ∈ M} the model space used at each nodal data point Y v .I.e.every node in the graph is associated with a model M v from a set common to all nodes, M. Thus a statistical model is associated with each node in the graph and hence each data point, The notation Θ v,Mv is used here to emphasis that since M v is the model order at node v, the parameter space will be dependent on it.For example, for compartmental models, the model order dictates the dimension of the parameter space.These collections of models can be used to generate a model space for the whole data set Y .Specifically, a set of candidate models for Y can be formulated as As mentioned before, the Potts model is a natural choice for the prior distribution over model order.For model selection, we adapt the hierarchical model above such that: 1. M = (M v ) v∈V is a Markov random field with a Potts distribution, with spatial dependence represented by E Y as before; 2. Each θ v , for all v ∈ V , is dependent on M v (e.g.model order determines parameter dimension in some cases, accordingly we denote the prior p(•|ξ, M v )) and known constant hyper-parameter ξ, importantly it is spatially independent (given the model order M v ); 3. Y v , the observed data, has likelihood f (•|θ v , M v ), where the model M v dictates which model is selected, for all v ∈ V .
Formally, the above can be summarised:

Discussion of Assumption 1
Assumption 1 may not hold in general: if there is spatial dependence between the generative model which describes each vertex then there may also be dependence between the parameters of those models.In some settings it can be viewed as an approximation which may be tolerable in order to allow inference under a model which is, at least, better than the assumption of full spatial independence.In other settings, this assumption may not be unrealistic.It is also noteworthy that, if the model order dictates the dimensions of the parameter space at that node, and the space contains dependent parameters, it may be difficult to incorporate spatial dependence at a the parameter level in a meaningful manner.
For example, if considering compartmental models: if two adjacent pixels were to contain different number of compartments, it may be difficult to say anything meaningful about the spatial dependence of the micro-parameters such as the transfer rates.Such considerations make it difficult to impose some more general Markov random field over the parameters.
Furthermore, incorporating spatial dependence at a model order level will typically take some precedence over spatial dependence at the parameter level.For instance, in the PET setting, generally we would like to know if two adjacent pixels have the same compartments (i.e.whether or not the localised area contains receptors) before we think about the transfer rates.
To summarise the problem of interest in general terms: the image Y is modelled node-wise using the density , with priors over the parameters θ and M .In particular, we propose the use of the Potts model as a prior over discrete parameters M to allow for tractable incorporation of spatial dependence in images.As mentioned above, θ will be treated as a latent variable -this is discussed in detail next.

Inference from the Proposed Model
Having specified this class of models, we now turn to the task of inference.Let p(M ) denote the Potts prior probability mass function over M ∈ M V .Then, the model posterior density, for Y , is denoted Here, f (y|M ) will be called the graph marginal likelihood of y; it encodes a generative description of the data under the proposed model.
The proposed method allows us to incorporate spatial dependence at the level of the discrete parameter M only.Subsequently, for computational tractability we impose Assumption 1: given the model order M v the random variable Y v at each node v ∈ V is independent.In other words, as per ( 5) and ( 6), we have that Here f (y v |M v ) denotes the likelihood at each node v, and will be henceforth called the node-wise likelihood.This assumption is important as it makes the graph likelihood tractable, and provides computational simplifications discussed later.Typically, as described above, a parametric model is used, thus f (y v |M v ) will in fact be a marginal likelihood or evidence (i.e. an integral, which is likely to be analytically intractable).
We finally have the probability mass function of primary computational interest Recall the mild assumption that hyper-parameter ξ is known and henceforth suppressed in notation.
Even in this general formulation, this mass function is not particularly attractive.There are some positives: It is straightforward to evaluate up to the intractable normalising constant; and Assumption 1 does seemingly lead to a simpler computation even at this stage, as only the node-wise marginal likelihoods f (y v |M v ) need to be computed rather than the high-dimensional integration of the full marginal likelihood f (y|M ).
However, even in this simpler setting, the graph marginal likelihood poses two computational difficulties: i) it is a product of |V | integrals and ii) these integrals in general are analytically intractable in most cases of interest.

The Node-Wise Pseudo-Marginal Algorithm
We now turn to extending the techniques discussed in Section 3, to build the novel computational method for use in the class of models constructed above.
The aim is to characterise the posterior density, π(M |y), of the model given in ( 7)-( 9).Suppose, for now, that a generic MH algorithm could be used in this setting.That is, to generate a MH Markov chain that targets the posterior density.This then gives the acceptance ratio: The NWPM Method Even with a tractable marginal likelihood, computing this acceptance ratio would be costly.Under Assumption 1 (given M the parameters θ = (θ v ) v∈V , of the model over the data y, are independent) the computational complexity of the problem is reduced considerably; now, only the node-wise marginal likelihoods need to be computed.However, every time a new configuration is proposed, up to |V | new integrals would then need to be computed and the acceptance probability is likely to be very small.Given that the number of possible Potts configurations is |M| |V | , growing geometrically with the number of nodes (i.e.pixels, for image data), the computational costs are not appealing.
Under Assumption 1 it is natural to employ a Metropolis-within-Gibbs approach to mitigate these difficulties and to obtain significant computational simplifications.Suppose that the MCMC sampler is currently at some configuration M = (M v ) v∈V .Under a single variable updating schedule: only the model order M v , at some node v ∈ V , will be proposed to be changed in a given step.We will refer to this as node-wise updating schedule, since we update the model order at a single node v only.Denote the proposed state to be M * v ; and M * [v] to be the vector equivalent to M but with M * v as the v-th coordinate.For simplicity, assume that the proposal distribution is symmetric, and the acceptance ratio is .
Thus, the acceptance probability of each Metropolis-within-Gibbs step requires the computation of only a single integral.
A complete sweep over the graph thus requires |V | such integrals but one could expect a reasonable proportion of these proposed moves to be accepted, in contrast to a global proposal which sought to update every node simultaneously.
Typically, even the node-wise marginal likelihood f (y v |M v ) will be difficult or impossible to evaluate analytically.A pseudo-marginal MH algorithm is a natural choice in such cases and is explored in the next section.

Graph Model Selection using Node-Wise Likelihood Marginal Estimates
The node-wise update schedule, gives rise to nested iterations.For clarity, term the outer iteration, indexed with i, of a full pass through the whole graph as the graphical iteration.The inner iteration over V will be termed the node-wise iteration.
Suppose that the MH Markov chain denoted, (M (i) ) n i=1 , is at graphical iteration i and that the state at node v is being proposed to be changed.Next, since the node-wise marginal likelihood is a normalising constant, for simplicity denote At graphical iteration i, at node v, for given model order proposal Associate with this generic estimator, the computation cost parameter N (for the SMC sampler this was the number of particles).As mentioned before this is a tuning parameter for the pseudo-marginal, thus also the algorithm presented below.
It is important to note the slight abuse of notation here - Following the node-wise updating schedule, whenever a new state ), is then .
We emphasise that the computation of r, precisely speaking, involves not just the standard MH proposal of a new configuration in M V , but also a proposal of the random variable Z v : it is a pseudo-marginal algorithm in the sense of Andrieu and Roberts (2009), with the added subtlety that the unbiased estimates of the marginal likelihood associated with every other node in the graph are also retained as a part of the extended state.The justification for using this acceptance ratio will be detailed in Section 4.2.2 The above choices and assumptions lead to the algorithm which we term the "Node-Wise Pseudo-Marginal" (NWPM) algorithm; presented in pseudo-code in Algorithm 2.
Algorithm 2 The Node-Wise Pseudo-Marginal Algorithm Given: (i) target density µ(M ) := π(M |y); (ii) unbiased node-wise marginal likelihood estimators With probability min{1, r} take: Within the context of this paper, the NWPM algorithm proposed here is a broadly applicable extension of the GIMH algorithm described above in Algorithm 1.
The output of this algorithm, the Monte Carlo sample (M (i) ) n i=1 , can be used to perform model selection in a principled manner.For example, this could be done using the node-wise marginal modal model order: i.e., at each node v ∈ V the model order which occurs the most in the (marginal) chain (M i=1 is selected.Additionally, samples of the latent variables θ are typically available as byproducts of the marginal likelihood estimates and can be used to perform parameter inference as demonstrated in the numerical studies below.
An advantage of this pseudo-marginal based approach in this setting is its flexibility.It can be used with essentially any unbiased estimator of marginal likelihood and hence allows existing technology from any application domain to be employed.The incorporation of spatial dependence via the modelling framework and algorithm developed here can come at the very little computational cost by approximating or making pragmatic adaptations to Algorithm 2. We employ such techniques in the simulation studies presented in Section 5 and discus in more detail in Section 4.4 below.
Finally, note that for convenience we assume that the tuning parameter N is fixed to the same value for all v ∈ V .This parameter typically determines the variance of the estimator Z v (M ) and additional flexibility could be obtained by allowing it to vary between vertices.Unsurprisingly, a larger N should give a better performance at the expense of greater computational cost.Indeed Andrieu and Roberts (2009), and more recently Andrieu and Vihola (2016), showed that when the dispersion of the marginal likelihood estimates are larger the mixing rate of any pseudo-marginal MH algorithm decreases.Of course, there is a higher computational cost when using larger N to reduce the estimate variance.The effects of this tuning parameter, and the trade off, on the mixing of the Markov chain in a standard pseudo-marginal algorithm have been studied by Sherlock et al. (2015); Doucet et al. (2015); Sherlock (2016).These works give theoretical justification and numerical studies for the recommendation that N should be chosen such that the variance of the log-likelihood estimator is close to one.

Marginal invariant distribution of NWPM
Here we establish the formal validity of a class of algorithms which includes the NWPM, allowing for more general moves than those described in Algorithm 2. We show that this class of algorithms leads to a Markov chain with an The NWPM Method invariant distribution coinciding with π, as given in (10).For further discussion on specifying proposals within this setting see Section 4.3.
Recalling Assumption 1, we may write the target distribution of the NWPM algorithm as Given M , denote the random vector of the unbiased node-wise marginal likelihoods estimators as " Z|M ∼ g(•|M ), i.e., Henceforth using the shorthand Z M v := Z v (M ) for all M ∈ M; we note that, due to independence, Finally, define the extended target density Here, we let f (y) denote the marginal probability of the observed data, integrating out unknown parameters and unknown model orders; Given observation Y = y, we treat it as a normalising constant.The NWPM algorithm has similar theoretical properties as GIMH.As such, the remainder of the argument of formal justification then follows the same as the GIMH case, see Andrieu and Roberts (2009).
Proposition 1 Let M denote a random variable in the graph model order space M V .Let " Z|M , g, µ and µ be defined as above.
For any U ⊂ V , let q U denote a Markov kernel on where q M U denotes a Markov kernel on M U and we slightly abuse density notation using the Dirac delta functions to indicate that variables associated with nodes outside U are unchanged (absolute continuity of the numerator and denominator of the Metropolis-Hastings ratio is ensured by the symmetry of this singular part of the kernel).
The standard Metropolis-Hasting acceptance probability with target distribution µ can be expressed as: and the marginal distribution of M under µ is µ.
Proof: First, consider the acceptance ratio of a MH (Metropolis-Hastings) algorithm with target distribution µ and proposal kernel q U : upon inserting the definition of q U one observes that the numerator is absolutely continuous with respect to the denominator and the singular elements simply impose that M −U = M * −U and " Z −U = " Z * −U , and we have: where the final factor is equal to one almost surely under the proposal distribution, and the result follows.
The marginal distribution of M follows by the essentially the same argument as in the standard pseudo-marginal context.
Clearly, the node-wise proposals described previously fit this framework with U = {u} being the single node being updated, although this framework would allow a somewhat broader class of proposals and blocked Metropolis-within-Gibbs type strategies to be explored.Standard arguments mean that a mixture or cycle of such kernels will also preserve µ as the invariant distribution.
One might anticipate that the theoretical properties of the NWPM algorithm might be further characterised by extending techniques from the standard pseudo-marginal setting such as Andrieu andVihola (2015, 2016); Andrieu et al. (2021) to this one, indeed it would be surprising were it to behave substantially differently, but doing this rigorously would involve some delicate technical work and is beyond the scope of this paper.
Next, we discuss some immediate extensions and innovations of this proposed methodology; resulting in further algorithms, approximations and future approaches.

Multiple Augmentation Pseudo-marginal Algorithms
In this section we will explore the potential of further augmentation of the state space.The specialised setting of the problem, allows approaches which to the authors' knowledge, have not been previously studied.The strategy developed here allows for a number of further extensions.Some simple approaches are discussed below -methods based on these approaches will also be evaluated in the numerical studies.
The augmentation of the state space used to justify node-wise pseudo-marginal algorithms can be further extended by adding an estimate of the marginal likelihood associated with every possible model at every node to the state space.
Although doing so may seem counter-intuitive and leads to a rather large state space, it allows a number of algorithmic innovations.In particular, by considering the following extended state space it will be possible to use a variety of standard MH moves in order to explore this space.Specifically, the further extended space allows us to decouple the updating of the marginal likelihood estimates from those of the model configurations themselves; the standard pseudo-marginal update would require that these two updates occur simultaneously.We exploit this decoupling when we go on to consider the possibility of updating the marginal likelihood estimates with lower frequency than the model configurations.
Let Z := ( Z v (M ) : v ∈ V, M ∈ M) be a vector of marginal likelihood estimator for each model order M ∈ M, at every node v ∈ V .Consider, next, a further extended form of the target density (12) : The extended joint density μ is constructed such that its marginal over M coincides exactly with the correct posterior distribution.Indeed, using essentially the same argument as in Proposition 1, we have: This further extended target distribution allows for some generalisations of the standard pseudo-marginal algorithm in the context of interest; it depends fundamentally on the fact that the variables associated with each node of the graph take values within a small finite set.It is possible to consider a variety of MH like moves applied to this extended target density.For simplicity we will consider only two types of proposal here: one which changes a single node's associated model and another which refreshes the likelihood estimates for a single node.
Denote the proposal of the random vector Z = ( Intuitively, Z * denotes a re-computation of all the components of Z.
First consider a proposal Q 1 u for the model associated with node u.
Same marginal estimates .
Here δ x,y is the Kronecker delta, taking value 1 if x = y and 0 otherwise and δ z (z * ) is, with the obvious abuse of notation, the singular measure concentrated at z evaluated over an infinitesimal neighbourhood of z * .In this proposal, only a change of M * u is proposed -in particular, Z and M −u do not change.Subsequently, the usual Metropolis-Hastings acceptance probability for this move is: .
For simplicity we have assumed that q 1 u (M u , M * u ) is independent of the remaining state variables, but this is not necessary for its correctness and proposals which depend upon their values could easily be implemented.Note that the augmentation of the state space ensures that Z M * u u exists and is available.
Similarly, consider a proposal Q 2 u which refreshes the augmenting variables associated with node u: for which the Metropolis-Hastings acceptance probability is simply by exploiting the same cancellation as in the standard pseudo-marginal setting.Clearly moves which update only some of the augmenting variables associated with node u could be justified in the same way.The augmenting variable Z M u for any M ∈ M can be re-estimated, not be just the normalising constant associated with the current state of M u .This flexibility is the result of the extension of the state space.
Standard arguments allow the combination of moves of these types, and others, within mixtures or cycles to provide an irreducible chain allowing considerable flexibility.In particular, it is no longer necessary to sample a marginal likelihood estimate for every proposed move in the state space.In the numerical studies below, we evaluate the performance of using combinations of such proposals.

Approximations of the NWPM Algorithm
We finish this section with a brief exploration of approximations of the exact pseudo-marginal algorithm described above with the aim of obtaining inference which is almost as good at a fraction of the computational cost by allowing for a small bias in those estimates.
Clearly, a significant portion of the computational load is used to compute the marginal likelihood; this is especially the case when using more sophisticated methods such as the SMC sampler.For the NWPM method, as presented in Algorithm 2, the pseudo-marginal MH Markov chain of length n requires n|V | marginal likelihood estimates.This quickly becomes infeasible with large data sets such as PET images, which may contain up to 10 6 voxel time series to be analysed.It is, therefore, worthwhile to consider strategies that reduce the number of marginal likelihood estimates that are required.
The multiple-augmentation approach to the NWPM algorithm decouples the estimation of marginal likelihoods from moves within the state space.One can envisage, for example, updating the marginal likelihood estimates only every κ iterations for some integer κ.Taking the limiting case as κ → ∞ and estimating the marginal likelihood once each model at each node suggests a cheap approximate scheme.Naturally, this approach leads to a Markov chain which is not irreducible on the extended space and which can no-longer be considered a pseudo-marginal algorithm.
Such an algorithm would simply sample each of the marginal likelihood of all the model orders only once, for each node, before starting the chain; then, follow Algorithm 2 otherwise, using only these initial single estimates.This amounts to making a stochastic approximation which will change the invariant distribution of the resulting Markov chain.
Doing leads to a (now marginal) MH Markov chain, that targets an approximation of the target density rather than the target density itself.As such, it would not be formally justified under the pseudo-marginal framework, and there would be a loss of accuracy with any subsequent results, since it is an approximation.Since we only do a single estimation of the marginal likelihoods we will refer to this method as the Node-wise Single-Estimation (NWSE) algorithm.
The biggest appeal of the NWSE is that, for a chain of the same length n, the sampler would only need to be used |V ||M| times, reducing the costs by n/|M| times.Any reduction in accuracy and performance could be adjusted for by using some of the saved residual computational resources towards reducing the variance of the sampler marginal likelihood estimates.
This approximation would allow at least preliminary spatial analyses to be conducted with little additional computational cost beyond that required for the associated mass-univariate analysis: if existing analysis has been done, and thus some estimates of the marginal likelihood have already been obtained then these can be readily used within the algorithm to incorporate spatial dependence.
The NWSE algorithm involves running an MCMC chain with the wrong invariant distribution: that in which the exact marginal likelihoods are replaced with the realisations obtained when they are sampled.It is of intereste, therefore, to establish how different such a distribution is from the true posterior distribution.The following elementary proposition, although it certainly does not provide a formal justification of the NWSE algorithm itself, demonstrates that in the context of small discrete spaces the error introduced by approximating marginal likelihoods can be controlled under reasonable conditions.
Proposition 2 For each model M ∈ M, let Z M := Z(M ) be a RV with expectation Z M := f (y|M ), corresponding to the associated marginal likelihood, and variance σ 2 M < σ 2 * < ∞.If the ( Z M ) M ∈M are mutually independent, then, letting M * = arg max{Z M } M ∈M which we assume to be unique the probability of selecting the correct model via maximisation of the marginal likelihood (equivalently, the posterior mode of the distribution over models given a uniform prior over M) is at least: by the Union bound.Applying Cantelli's inequality, noting that Z M * − Z L > 0: Combining these expressions yields: This bound can be made arbitrarily close to 1 by choosing estimators with sufficiently small variance.As the variance of the normalising constant estimates needs to be assessed to allow pseudo-marginal algorithms to be tuned, it is reasonable to suppose that this information could be obtained in settings in which this type of method is used and hence that a reasonable degree of confidence can be obtained that the use of this approximation does not substantially influence the resulting inference.Naturally, this result suggests that if estimating each marginal likelihood once rather than for every algorithmic step as in a pseudo-marginal algorithm, a relatively small variance might be required to obtain good performance.
In the numerical study, presented below, we investigate, evaluate and compare the effects when using the NWSE algorithm in different settings.

Applications
In this section, the NWPM algorithm and some variants are applied to a simple toy model, in which the true marginal likelihood is known, and to a compartmental model for PET data with both simulated and real data.
The software implementation of these methods can be found in the R package bayespetr, available at: https: //github.com/dt448/bayespetr.The proposed algorithms are somewhat computationally intensive and here we focus on rather short Markov chains, a situation which we expect to be most widely useful; in Appendix B we establish that running these chains for much longer does not materially change the resulting inference (as we have found to be the case in all settings which we have explored).
The 20 × 20 image displayed in Figure 2 will be used to generate the ground truth Potts configuration for all simulateddata experiments.It was adapted from Bezener et al. (2018) and includes a range of spatial structures of varying complexity.The image is split into four regions, labelled R 0 , . . ., R 3 ; all pixels within a given region have a common model order; the details vary between experiments and are given below.
Whenever the NWPM algorithms are used in the numerical studies below, M was initialised using the Gibbs a sampler targeting the prior (Potts) distribution.However, when analysing measured PET data in Section 5.3, the algorithms were initialised from the output of the spatially independent SMC sampler method.
In the empirical studies below, we investigated a simple proposal strategy based on the multiple augmentation space, as discussed in Section 4.3.
In brief, the proposal kernel Q 1 , together with Q 2 at every κ ∈ N graphical iteration, is used.Essentially, upon initialising the chain, only the change of the model at the node is proposed.After every κ−th graphical iterations, changes in the auxiliary variable or marginal likelihood estimates for every model order at every node is proposed to change.In this section, we will refer to this variant as the NWMA algorithm, with tuning parameter κ.

SMC Sampler:
The algorithms considered require unbiased estimates of f (y The SMC sampler can be used to compute these node-wise marginal likelihoods by targeting the parameter posterior distribution π(θ v |y v , M v ) using an annealing schedule.
More precisely, for each model order M ∈ M, let the sequence of target distribution of the SMC sampler be {π t } t≥1 ; Here, we define ) , where the total number of intermediate distributions T and annealing scheme, α : [0, 1] → [0, 1], may be different for each model order (Zhou et al., 2016, Algorithm 2).For example, the annealing scheme could be a simple fixed schedule such as α(t/T ) : [0, 1] → (t/T ) 5 , which following Zhou et al. (2016), we refer to as the "Prior 5" annealing scheme.Alternatively, α t can be determined adaptively based on metrics of the distance between intermediate distributions, such as the ESS (Jasra et al., 2011;Kong and Liu, 1994) or CESS (Zhou et al., 2016).
For the numerical studies below, we use the SMC sampler.Pilot studies showed that Prior 5 and the CESS-adaptive annealing scheme produced similar results in these examples, thus the Prior 5 scheme is used for simplicity.

Toy Model Simulation Studies
Toy Model: A simple toy model in which both the prior and model likelihood are normal, at every node, is used to validate the proposed methodology.The data-set comprises a 2-dimensional digital image, each pixel having a scalar intensity.The graph G = (V, E) used to represent the spatial structure is a finite square lattice, i.e., it has vertex set In the notation of Section 4.1.2,given an image Y = (y v ∈ R : v ∈ V ), the model studied here is: where N (µ, σ 2 ) denotes the normal distribution of mean µ and variance σ 2 > 0.
At each pixel, v, we model the data point Y v as a normal random variable with mean µ v and variance σ 2 , which is assumed known and is common to all pixels.The prior over the latent variable µ v is normal with mean µ (Mv) 0 determined entirely by the model order M v and fixed variance σ 2 0 .As before, the model order configuration will be the random variable M with a Potts prior distribution.
The node-wise marginal likelihood can be straightforwardly evaluated: , for all M ∈ M.This is trivial to compute since all these terms are known.
The spatial configuration of model order used to simulate the data was specified to represent simple spatial structures; the ground truth value of M was fixed to be the Potts configuration shown in Figure 2. Finally, we specified σ 2 0 = 5 2 and σ 2 = 1 2 .

Simulation Study 1: The Coupling Constant
As discussed in Section 2.3, typically it is difficult to infer the coupling constant J of the Potts distribution (Møller et al., 2006;Moores et al., 2020).In this work, we make the assumption that J is known when using the NWPM algorithmthis is not too unreasonable as J is a parameter of the prior distribution.In this vein, here we evaluate the performance of the algorithm for various values of J.
The following design was used: An image data set of 20 × 20 = 400 pixels(nodes) from the toy model and model order ground truth as dictated by the Toy Model Ground Truth above, was generated.This simulated image was then analysed using the NWPM algorithm, for varying values of the coupling constant J. Coupling constants ranging from J = 0 to Figure 3: Empirical average, and standard deviations (s.d.), of the percentages(%) of the total (400 or 10,000) pixels where the correct model order was selected, for different values of the coupling constant J. Model selection was done on simulated toy model data using NWPM with an SMC estimator (with N = 50, and T = 80) for the marginal likelihood.The configuration selected for estimation by looking at the percentage of time each node spent in each model order state.Both line plot shows the empirical average of the correct percentage of pixels.The lighter and darker region show the 2× s.d. and 1× s.d.error bands, respectively.J = 5 were investigated.Based on pilot studies, an SMC sampler with N = 50 particles and T = 80 distributions was used within the NWPM algorithms.The pseudo-marginal Markov chain was ran for n = 100 iterations.As mentioned above, the pseudo-marginal MH Markov chain was initialised using a Gibbs sampler, with the respective coupling constant.This was done for 50 replicates of the NWPM Markov chain.
Model selection was carried out by selecting the modal model order marginally at each node v ∈ V .In other words, the modal state in which each M (i) v , for i = 1, . . ., 100 was selected for each pixel v ∈ V .The empirical averages, over the 50 replicates was calculated.The mean percentages of the nodes for which the correct model order was selected for the different variations of the NWPM algorithm is shown in Figure 3a.
It is evident from looking at Figure 3a the performance of the algorithm peaks for coupling constant J = 0.4.This is smaller than the critical threshold J critical , as discussed in Section 2.3.This is noteworthy, as there is typically a sharp decrease in speed of mixing of the Markov chain above this critical value.In fact, the graph shows that the variance of the percentage of correctly selected model orders increases with the coupling constant; Importantly, there is a jump in the rate of increase between J = 1 and J = 1.5.
Given that the size of measured PET data will be considerably bigger, this experiment was repeated for a larger Toy Model image.Additionally, the model order space for compartmental models may be larger; For instance, in this paper we consider up to M = 3−compartment models.Subsequently, a scaled version of the 20×20 image, with dimensions 100×100 = 10,00 pixels and model space M = {C, D, E}.More specifically, the region corresponding to R 0 had mean parameter µ (C) 0 = +7; Similarly, region R 1 had µ (D) 0 = 0 and regions R 2 and R 3 had µ (E) 0 = −7.The parameters σ 2 0 and σ were kept as above.The generated data was the analysed using the same the methods and design.The results are shown in Figure 3b.
Notably, we first see that the peak of the graph is flatter in comparison to Figure 3a.Secondly, the performance of algorithm is better overall; For instance, the best performance is almost 100% compared to roughly 94% previously.Similarly, when looking at the performance corresponding to J = 0.4, it increases from 94% to 97%.This suggest that, in this context, a larger image size and model order space did not seem to impact the performance negatively.Similarly, there is evidence that a wider range of coupling constant values could be used in the setting with bigger images and larger model order space, to attain "optimal" model selection performance.A selection of model order outputs, for varying values of J, are given in Figure 17.
When looking at these results, as presented in Figure 3, an important question to consider here is how much of the performance difference is due to poor approximation of the posterior by the Markov chain, and how much to the posterior itself.It is reasonable to argue that below the critical value, reasonably good mixing is obtained and this probably reflects a good approximation of the posterior as indicated by the low variability.In contrast, above this critical value the results are largely dominated by the poor mixing of the Markov chain over this rather short period.We also note that the graph shows good performance for even this relatively short chain, for the appropriate coupling constant values.
Based on evaluating the performance for a range of coupling constant and due to the similar qualitative features of the different settings, in the interest of simplicity, for all the following studies we will fix J = 0.4.This lower value means that the risk of being too close to the critical value J critical is avoided.

Simulation Study 2: Sample Size vs Chain length
Next, the computational trade-off between the number of particles used in the SMC sampler and the length of the pseudo-marginal MH Markov Chain was investigated.The same design as the study in Section 5.1.1 was used, with the exception of the Markov Chain length varying from n = 50 to n = 200, and the number of particles varying between N = 50 and N = 200.Here, we investigate the 20×20=400 pixel image.
Additionally, spatially independent model selection, using an SMC sampler to estimate the model evidence was also investigated.The NWSE approximation was also used to analyse the image.Similarly, the NWMA algorithm, as described above, with updates to the marginal likelihood estimates every κ = 10 graphical iterations was used.For this variant, an SMC sampler with N = 200, T = 500 was used.For these two, less resource-intensive, Markov chain methods (NWMA and NWSE), n = 200 graphical iterations were used.
Model selection was carried out in the same manner as described above (i.e. the modal model), and the percentage of correctly selected nodes was used as the metric of performance.This was done for 100 replicates, for each method used.Results are shown in Figure 4.These results suggest that, at least in this simple setting, the algorithm is not particularly sensitive to the allocation of computational resources.Importantly, we can see that the selected graph model/configuration stabilises fairly quickly for all the combination of tuning parameters.In particular, looking at the longer n = 200 chain in Figure 4a, it is evident that the convergence is to around 92% of nodes having the correct model order selected when using the NWPM algorithm.This is considerably higher than the performance of the spatially independent SMC sampler, which was 79%; confirming that exploiting spatial structure can significantly improve estimation in this type of model.We also observe that the NWSE approximation method and the simple NWMA algorithm achieved almost identical results to the full NWPM algorithm: suggesting that, in this simple setting, there is enough signal to allow for the use of the approximation to achieve very similar results, with very little additional computational cost.Intermediate chain lengths not shown in the graph, i.e. n = 75 and N = 134, also produced similar results.
Figure 4b shows the average computational runtime for each method.The runtimes were computed by considering the average CPU time per iteration: That is, the total run time for each trail was recorded and then averaged over the total graphical iteration.We can draw similar conclusions as above, with the exception that the NWPM method with N = 200 particles seemed to take a longer time (in minutes) to reach stability.The computational runtime for the NWMA method, for varying values of κ, was also studied: Results are shown in Figure 15.As expected, the total runtime increase for smaller values of κ; However, a higher refresh rate seems to slightly slow down (in minutes) how quickly the method attains its stable state.
Finally, the long term behaviour for the NWSE and NWMA MCMC chains was also explored and results (i.e. that much longer chains lead to very similar conclusions) are shown in Appendix B, Figure 13.Additionally, Figure 11 shows an enhanced version of Figure 4 -allowing for better discernment of the small differences in performance.

Simulated PET Data Studies
Before we investigate the performance of the proposed algorithms on measured PET data, we verify and investigate tuning factors when applied to simulated PET data.
Synthetic data was simulated using the same plasma input function and integration periods as in the measured PET data set studied in Section 5.3 and normally-distributed noise, with zero mean, was added.The variance was such that it was proportional to the true time activities divided by the length of the frame duration; and scaled such that the highest variance in the sequence is equal to the noise level.A noise level of 0.5, which lies within the representative range (between 0.01 and 1.28) of real data (Jiang et al., 2009), was used when simulating the 2-D PET image.
We employed vague uniform priors over φ and ϑ, which were constrained to lie within [10 −5 , 10 −1 ] and [10 −4 , 10 −1 ], respectively.These ranges are based on pilot and previous studies, and ensure that the parameters are physiologically meaningful (Cunningham and Jones, 1993).
The precision parameter of the measure noise is denoted λ = 1/σ 2 .A gamma distribution, with both parameters equal to 10 −3 , was used as the prior for λ -a proper approximation to the improper Jeffrey's prior.
Details of the relevant distributions can be found in Appendix A.
PET Simulation Ground Truth As before, the ground truth of the model order configuration were based on spatial structure pattern of Ground Truth Configuration, shown in Figure 2. In this PET simulation study, a 20 × 20 Potts configuration with model order space M = {1, 2, 3}, representing the number of compartments was generated.Region R 0 corresponds to model order 2; region R 2 and R 3 to model order 1 and region R 1 to model order 3.

Simulation Study 1: Algorithm Comparison
Pilot studies were initially carried out to evaluate the model selection performance and estimator variance of the SMC sampler.For the variance of the estimates, on agreement with numerical studies by Zhou et al. (2016): the numerical study showed that, there was a decrease linearly for higher values of N and similarly in T .However, beyond a certain threshold, the reduction were no longer useful.Similar results were seen in pilot studies for model selection performance, for 1-D data simulated from just the 2−compartment model with noise level 0.5, performance could not be increased beyond 65% regardless of any increase in tuning parameter values.
Next, the simulated 2-D PET image was analysed using the NWPM algorithm, and the effect of the pseudo-marginal MH chain was investigated.Based on the above pilot study, an SMC sampler with the Prior 5 annealing scheme, T = 400 intermediate distributions, was used as the estimator of the normalising constants at all pixels for all models.These values, and the values of the number of particles N in the design below, were chosen as they allow for the variance of the likelihood estimates to be within a suitable range to allow for optimal scaling (Doucet et al., 2015).
Chain lengths of n = 50, n = 75, n = 100 and n = 200 were investigated, with particles N = 200, N = 134, N = 100 and N = 50, respectively.These, relatively smaller, graphical iteration lengths were chosen due to the computational overhead of this algorithm in this resource intensive setting.Note that, since model selection here uses the marginal modal state from a space of just 3 model orders -these shorter chains can still produce reliable results.
Additionally, the image was also analysed using spatially independent SMC sampler and the NWSE approximation; using an SMC sampler with N = 400 particles and T = 600.The NWMA algorithm, with SMC sampler using N = 200 and T = 600, at every κ = 50 iterations, was also evaluated.Given their significantly reduced computational costs, for the two pseudo-marginal methods, n = 500 iterations were used.
Following Algorithm 2, all the chains where initialised using the Gibbs sampler to target the prior distribution.30 independent replicates of each algorithm were used to obtain the results shown in Figure 5 and Table 1 (performance for the intermediate ranges of tuning parameters of the NWPM algorithm is shown in Appendix B, Figure 16).
Looking at Figure 5a, the results show notable improvements in model selection performance when using the different variants of the NWPM algorithms compared to spatially independent model selection using the SMC sampler.Analysing the simulated 2-D PET image using the spatially independent method gave 61.4% pixels correctly selected, compared to the NWPM algorithms all giving around 76%, even for as short as n = 50 graphical iterations.On average, this is a 23.4% increase in the number of pixels selected correctly when using the proposed NWPM approaches compared to the spatially independent SMC sampler method.On large data-sets, like PET images, this is a considerable number of pixels.Furthermore, since this is spatially dependent model selection, the improvement in correctly selected pixels when compared to SMC independent will be towards revealing the underlying spatial structures.These results suggest that incorporating spatial dependence does results in better model selection performance.
Figure 5 also shows that the performance of the NWSE approximation and the NWMA variant both provide similar performance to the full NWPM algorithm.Pragmatically, this is very promising as it suggests that even with an approximation current state of the art methods can be outperformed with very little additional costs.Alternatively, rather than using an approximation, algorithms based on the multiple augmented state space with significantly reduced computational overhead can be used.For example, the simple NWMA algorithm studied here.Similar to the toy model, there does not seem to be any notable difference when using higher chain length and lower particles.Additionally, there is only a small increase in performance when using more a computationally intensive algorithm, on average.That is, the full NWPM performs slightly better than the NWMA variant; which in turn also has a minute increase in performance compared to the NWSE approximation.
Additionally, for the chains where higher iterations were used, we can see that this increase in results was maintained.To verify long term behaviour of the algorithms: long-run MCMC chains, of the less expensive NWSE and NWMA algorithms, are shown in Appendix B, Figure 14.As before, an enhanced version of Figure 5a is shown in Figure 12.
As before, Figure 5b shows the average computational runtime for each method.The runtimes were computed in the same manner as described in Section 5.1.2.Similarly, as seen in the aforementioned section, the NWPM method with n = 200 seems to be the slowest in regards to mixing (in minutes).
Next, Table 1 shows the average mean squared error (MSE) of the V D estimates for each variant of computational method.The quantity was calculated by taking the average MSE over all the pixels in the 2-D PET image, and then over all the replicates (the empirical variance of these quantities are shown in the parenthesis).Although the non-spatial SMC sampler method has the highest empirical MSE, the NWPM algorithms, including all chain lengths and the SE variant, produce a noticeably smaller MSE compared to the spatially independent SMC method, and there is a small improvement when using the full pseudo-marginal NWPM, compared to its approximation, these differences all fall within the range that could be plausibly explained by random fluctuations and it is difficult to draw any meaningful conclusion.

Measured [ 11 C]diprenorphine PET image Study
The simulation studies above have verified the effectiveness of the proposed methodology and allowed for the identification of tuning parameters which lead to promising results in settings close to that of interest.We now turn our attention to real data sets.
Analysis of the tracer kinetics of dynamic PET images, that use the [ 11 C]-diprenorphine tracer, are generally aimed towards quantification of opioid receptor concentrations in the brain.Such studies involve investigation of neurological disease which often involve change in brain receptor density.
Investigation of similar PET data sets, using different statistical approaches, have been previously analysed quantitatively by Gunn et al. (2002), Peng et al. (2008) and Jiang et al. (2009); However, these works focused on parameter estimation rather than model selection.Zhou et al. (2013Zhou et al. ( , 2016) ) are more recent works that investigated model selection; Albeit, in these works, spatial independence is assumed.

PET Data
A test-retest pair of measured dynamic PET data of the concentration of the tracer [ 11 C]-diprenorphine in a normal subject, for which an arterial input functions were available, was analysed using the proposed and other methods.This data was initially measured as part of the dynamic PET study in Hammers et al. (2007).In that study, the effects of using lower concentration of tracer was investigated.In addition, the impact of using different slow frequency cutoffs within the analysis of the data was also studied.
PET data acquisition The subject underwent 95 min dynamic [ voxels.This gives a total of 1, 556, 480 separate times series respectively to be analysed.In this study, we will analyse a cross-section of each of the sagittal, coronal and transversal planes of the brain -the total number of times series studied here is, therefore, roughly 11,000.
Details of further PET data pre-processing can be found in Hammers et al. (2007, Material and methods).

Decay Correction
The PET data is not decay corrected and this should be accounted for within the compartmental model.Gunn et al. (2002) (similarly, Cunningham andJones (1993) for Spectral Analysis), suggest that the slow frequency cutoff be close to the decay constant of the tracer used.In other words, given that the half life of [ 11 C] is 20.4min, its decay constant is 0.0005663s −1 ; Thus, the slow frequency cutoff, denoted θ min , is fixed close to this constant.Hammers et al. (2007) showed that actual parametric values depended heavily on the cutoff slow frequencies( between 0.0008 s −1 and 0.00063 s −1 ).Importantly, the prior ranges of θ 1:M should be based upon θ min , and subsequently this should be taken into consideration when computing the volume of distribution V D .Following Zhou et al. (2016), in particular see Zhou (2014), we will use the cutoff 0.0007s −1 .

Data analysis and performance evaluation
Normally-distributed errors have been found to poorly model data of this sort (Zhou et al., 2013).Instead we will use the t−distribution to model the additive error variable .The same gamma prior as used for λ in the simulation studies above, will be used for the scale parameter, τ .A uniform distribution over the interval [0, 0.5) will be used as the prior for 1/ν, allowing the likelihood to vary from having a very heavy tail to being arbitrarily close to normality.
The measured data was analysed using the following methods: Firstly, SMC sampler estimates of the marginal likelihood of each pixel for the three models were computed.This strategy assumes spatial independence, so, as before, we will refer to this method as the independent SMC method.The Prior 5 annealing scheme with N = 300 particles and T = 500 distributions was used, these values were based on the studies above and numerical results reported by Zhou et al. (2016).Three cross-sectional slices (coronal, transverse and sagittal planes) of the 3-D PET image were analysed, the model order and volume of distribution parametric images are shown in Figure 6.
Next, the NWPM algorithm was used to analyse the same three cross-sections.An SMC sampler with N = 300 particles and T = 500 distributions was used as the marginal likelihood estimator.Here, the chain was initialised from the output of the SMC spatially independent model selection above (rather than from the prior distribution, as in the simulation studies above).A pseudo-marginal MH chain of length n = 75 was generated.This chain length was based on viability of computation cost for this relatively large data set, and the simulation studies above.In addition, the sagittal cross-sectional image was analysed using a longer n = 200 iterations, as discussed below.The results are shown in Figure 7.
Finally, Figure 8 and 9 show the model configuration output from analysing the measured PET data using the NWSE approximation and the NWMA variant, respectively.In this case, an SMC sampler with N = 400 particles and T = 600 distributions was used to make the marginal likelihood estimates.The pseudo-marginal chain was initialised using the output of these estimates, and these initial estimates where re-used within the NWSE algorithm for n = 500 iterations.
For the NWMA method, the chain was initialised in the exact same manner, but the marginal estimates where updated every κ = 100, for a total of n = 500 iterations.
For further comparison, the sagittal cross-section was also analysed using NLS (non-linear least squares; Zhou et al. (2013)).This is shown in Figure 10.The spatially independent SMC sample model selection and NWPM method(with a longer chain length n = 200) model selection are included to show the difference in the final model order parametric image output.
The volume of distribution images, specifically Figure 6 and Figure 7, show very similar performance for all methods when considering inference of V D .However, Figure 10 shows clear difference in the model order image when using NWPM compared to spatially independent SMC and NLS.It is evident that analysis using NWPM reveals more spatial structure: Albeit, using the SMC independent approach does reveals some of the underlying spatial information when compared to NLS; NWPM has a de-noising effect and improves the clarity of the roughly formed clusters seen in the spatially independent SMC method.In particular, the M = 1 and M = 2 compartment models are selected to model most of the structural clusters.In contrast to this there is a decrease in the number of M = 3 compartment models selected when compared to the non-spatial SMC method.
Note also that, by looking at Figure 10 we see that the NWPM algorithm converges to a noisy version of the final output image even at low iterations as n = 25.Importantly, there seem to be only small changes at the higher iterations, suggesting that the chain reaches a stable output relatively quickly.Similar behaviour is seen in the simulation studies above, albeit in a smaller, simpler settings.
Furthermore, similar model selection output is seen when using the SE estimate approximation and the multiple augmentation variant, NWMA, of the NWPM algorithm.Comparing both Figure 8 and 9 with Figure 7, we see that the mode order outputs are not identical, but qualitatively show many similarities.In particular, we see that the there is a decrease in the total number of 3-compartment models selected when using the SE approximation or the NWMA variant, compared to the full NWPM.A similar effect is seen when comparing the independent SMC model selection with the NWPM.Notably, the NWSE and NWMA model selection outputs are almost identical, differing only at a very small number of pixels.

Conclusion
In this work we have illustrated a novel computational method for effectively incorporating spatial dependence when performing model selection at each location within a graph.This approach extends the pseudo-marginal method, in a number of directions, within the context that it has been developed and evaluated in.By exploiting the structure of the problem at hand, the developed methodology allows for considerable flexibility in the updating of state variables and marginal likelihood estimates relative to generic pseudo-marginal algorithms.The empirical study suggests that this approach can yield better inference than methods which impose assumptions of full spatial independence and, in particular, can reveal clearer spatial structure in the image of underlying model orders.As with pseudo-marginal algorithms, essentially any unbiased marginal likelihood estimators used in non-spatial analysis can be used within this algorithm.
This methodology was motivated by considering the problem of model selection in PET images.In such a setting both the amount of data and the complexity of the models lead to substantial computational requirements.The images of model order provided by this approach supplement the information provided by volume of distribution and similar macroparameters in PET images.For instance, Figure 10 (part d.) could be used to make the interpretation that the signal from some parts of the PET image require a higher model order compared to other parts .Specifically, this is apparent when comparing the posterior region of the brain image (where 2-compartmental models were selected) to the anterior region of the brain (where 1-compartmental models were selected).
In order to reduce the computational cost, a multiple augmentation variant of the pseudo-marginal algorithm is introduced and a more extreme approximate form in which each marginal likelihood is estimated only once.The performance of these methods in the examples explored here is very encouraging: it suggests that where non-spatial The NWPM Method analysis has been performed, spatial dependence can be incorporated with the existing estimates very easily and with little computational costs.Doing so carefully can result in similar performance to the full NWPM algorithm.
In summary, the NWPM algorithm, together with its approximations and extensions, enables a intuitive, natural approach to incorporating spatial dependence in realistic and complex settings; such as analysis of PET images.In particular, the methods extensions and approximations means that the subsequent performance and results of incorporating spatial information is very accessible and readily applicable.

Figure 1 :
Figure 1: Plate diagram of the proposed generic model.
Bayesian model selection of spatial data Y can be thought of as inference of the model order parameter M in the space M V .Each realisation of M is called a configuration, note that there are |M| |V | candidate models for Y .

Figure 2 :
Figure 2: Ground Truth Configuration: spatial structure used to produce ground truth Potts configuration for the simulated data studies.A 20 × 20 grid, with four homogeneous regions.

Toy
Model Ground Truth: Ground Truth Configuration, as shown in Figure 2 was used to generate a ground truth model order image.For this toy model, this is a 20 × 20 Potts configuration with model order space M = {A, B}.Region R 0 was fixed to model order M v = A and the remaining three regions R 1 , R 2 and R 3 where fixed to model order M v = B.At each node v ∈ V , the mean parameter µ (Mv) 0 will depend on M v .In other words, the lighter pixels have nodes with hyper-parameter µ (A) 0 = +5 and the darker pixels have nodes with hyper-parameter µ Model selection performance for 100×100=10,000 pixel image.
Correctly Selected Node Models for varying MCMC chain lengths and number of SMC particles (a) Mean percentages at each graphical (MCMC) iteration.average CPU runtime using varying MCMC chain lengths and number of SMC particles (b) Mean computational runtime (averaged over each graphical iteration).

Figure 4 :
Figure 4: Average percentages (%) of the whole image (20 × 20 = 400 pixels, Toy Model) where the correct model order was selected at each iteration (shown in part (a)) of the MCMC chain, using NWPM, NWSE and NWMA for varying Markov Chain length, n, and number of particles, N, in the SMC sampler.The dashed line shows the average percentage when using spatially independent SMC model selection.The cumulative average over the graphical iterations i = 1, . . ., n is shown.Average computational runtimes are shown is shown part (b).
Nodes Correctly Selected for Simulated PET data(noise level = 0.5) using different methods.(a)Mean percentages at each graphical (MCMC) iteration.average computational time using varying MCMC chain lengths and number of SMC particles for Simulated PET data(noise level = 0.5) (b) Mean computational runtime (averaged over each graphical iteration).

Figure 5 :
Figure5: Average percentages(%) of the whole image (20 × 20 = 400 pixels, simulated PET data) where the correct model order was selected at each iteration of the pseudo-marginal MH chain, using: NWPM for varying Markov Chain length, n, and number of particles, N, in the SMC sampler; NWSE with N = 400 particles and T = 600 distributions and NWMA using SMC sampler with N = 200 and T = 600 at every κ = 50.The dashed line shows the average percentage when using spatially independent SMC(N = 400, T = 600) model selection.

Figure 12 :
Figure 12: Average percentages(%) of the whole image (20 × 20 = 400 pixels) where the correct model order was selected at each iteration of the pseudo-marginal MH chain.MCMC traces, of each method, for graphical iterations n = 10 and on wards are shown.
Figure 14: MCMC traces of the percentage of correctly selected nodes, for simulated 2-D PET image, using the NWMA (refreshing marginal estimates at every κ = 1000 graphical iterations) and NWSE methods.For both methods, the SMC sampler with N = 200, T = 400 was used.

Table 1 :
Average RMSE (and s.e.) of V D estimates for simulated PET data.
(Kinahan and Rogers, 1989)eline scan.[ 11C]-re norphine is a tracer, of the non-selective antagonist, which binds to neural opioid receptors.The PET scans were acquired in 3D mode on a Siemens/CTI ECAT EXACT3D PET camera.After image reconstruction, the spatial resolution is approximately 5mm.Technical details including the data correction and normalisation methods of the ECAT EXACT3D PET camera can be found inSpinks et al. (2000).Thirty seconds after the start of the scan, [ 11 C]-diprenorphine was injected and flushed through the cannula as a smooth bolus over a total of 30s.The subject was injected with 185 MBq (Mega Becquerel) of [ 11 C]-diprenorphine.The PET data was reconstructed using the re-projection algorithm(Kinahan and Rogers, 1989)with ramp and Colsher filters cutoff at the Nyquist frequency.Reconstructed voxel size were 2.