2.1. Extending the basic model–mixture formulation

We now extend the model by proposing a Poisson mixture model formulation. We assume that D _i, given λ _j, has a Poisson distribution with mean λ _jE _i. We take a non-parametric Bayesian approach to modelling this mixture distribution using a Dirichlet process prior, and use RJMCMC to estimate the number of components in the mixture. In this case, the physical interpretation of the model is that the heterogeneity is assumed to be drawn from one of k possible components, in proportions w ₁, … , w _k. The method we describe is essentially a classification problem where we assume that each observed D _i comes from only one of k components, where each component has a Poisson distribution. Thus,

$$D_{i} \mid\lambda _{j} \sim{\rm Poisson}\,(E_{i} \lambda _{j} )\quad j\,=\,{\rm 1},\ldots,k;\,i\,=\,{\rm 1},\ldots,n$$

More general forms of the mixture Poisson model with covariates are discussed in Green & Richardson (Reference Haastrup2002). Mixture models for grouped claim numbers are considered by Tremblay (Reference Tremblay1992) and Walhin & Paris (Reference Walhin and Paris1999, Reference Walhin and Paris2000). Dellaportas et al. (Reference Dellaportas, Forster and Ntzoufras1997) considers count data in finance using split/merge moves, whereas Viallefont et al. (Reference Viallefont, Richardson and Green2002) provide a more general discussion of mixtures of Poisson distributions using both split/merge moves and birth/death moves. Other methods for determining the number of components in a mixture are discussed by McLachlan & Peel (Reference Ntzoufras and Dellaportas2000), Phillips & Smith (Reference Phillips and Smith1996), Carlin & Chib (Reference Carlin and Chib1995), and Stephens (Reference Viallefont, Richardson and Green2000), who use Markov chains to model jointly the number of components and component values. The advantage of the Bayesian formulation is that we can place posterior probabilities on the order of the model.

2.2. The likelihood function

Throughout our discussion, n will denote the number of data points and k the number of components in the mixture formulation. For a finite mixture model, the observed likelihood function is

(1) $$L\,(D^{n} \mid\lambda ,\,w,\,E^{n} )\,=\,\prod\limits_{i\,=\,1}^n {\mathop{\sum}\limits_{j\,=\,1}^k {w_{j} f_{j} \,(D_{i} \mid\lambda _{j} ,\,E_{i} )} } $$

where the weights are non-negative and $\mathop{\sum}\nolimits_{j\,=\,1}^k {w_{j} \,=\,1} $ . Even for moderate values of n and k, this takes a long time to evaluate as there are k ⁿ terms when the inner sums are expanded (Casella et al., Reference Dellaportas and Karlis2000). Another form of the likelihood function will be derived shortly. Classical estimation procedures for mixture models are described by Titterington et al. (Reference Walhin and Paris1990) and McLachlan & Peel (Reference Ntzoufras and Dellaportas2000).

Let z _i be categorical random variables taking values in {1, … ,k} with probabilities w ₁, … ,w _k, respectively, so that

(2) $$p\,(z_{i} \,=\,j\mid w)\,=\,w_{j} $$

Let f _j (·) denote a Poisson density with parameter λ _j. Suppose that the conditional distribution of D _i, given z _i=j, is Poisson (λ _j), j=1, … , k. Then the unconditional density of D _i is given by

$$\eqalignno{ f\,(D_{i} )\,= & \mathop{\sum}\limits_{j\,=\,1}^k {f_{i} \,(D_{i} \mid z_{i} \,=\,j,\,\lambda ,\,E^{n} )\,p\,(z_{i} \,=\,j)} \cr = & \mathop{\sum}\limits_{j\,=\,1}^k {w_{j} f_{i} \,(D_{i} \mid\lambda _{j} E_{i} )} $$

With each pair (D _i, E _i), we associate a latent variable z _i, which is an indicator variable that indicates which component of the mixture is associated with (D _i, E _i). We have z _i=j if the ith data point (D _i, E _i), comes from the jth component of the mixture. Thus, for each i, we have

$$z_{i} \mid w\sim{\scr M}({\rm 1};\,w_{{\rm 1}} ,\ldots,\,w_{k} )$$

and

$$D_{i} \mid z_{i} \sim{\cal P}(\lambda _{{zi}} E_{i} )$$

By incorporating the indicator variables z _i, the complete data likelihood is then

(3) $$\eqalignno{ L(D^{n} \mid z,\,\lambda ,\,E^{n} )\,= & \prod\limits_{i\,=\,1}^n {f\,(D_{i} \mid\lambda _{{z_{i} }} ,\,E_{i} )} \cr = & \prod\limits_{j\,=\,1}^k {\prod\limits_{\{ i:z_{i} \,=\,j\} } {f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )} } $$

At times, especially for the fixed k case described below, it is more convenient to work with (3) as it involves multiplications only, rather than additions and multiplications, as in (1). Note that the inclusion of the categorical variables (z _i) does not add to the complexity of the model. Instead, it results in a simplification of the likelihood function from being a product of sums as in (1) to a product only, as in (3).

The convenience of using the missing data formulation is that the posterior conditional distribution of the model parameters would be standard distributions. Moreover, the augmented variables z _i allow us to see the component of the mixture to which the data points are assigned.

2.3. Gibbs updates for fixed k

We consider a mixture of Poissons, where conditional on there being k components in the mixture:

$$D_{i} \sim\mathop{\sum}\limits_{j\,=\,1}^k {w_{j} f\,( \cdot \!\mid\lambda _{j} ,\,E_{i} )} $$

The weights w _j, sum to 1, and are non-negative. That is,

(4) $$\mathop{\sum}\limits_{j\,=\,1}^k {w_{j} \,=\,1\;{\rm and}\;w_{j} \,\geq \,0} $$

f(D _i ∣ λ, z _i=j, E _i) ~ Poisson (λ _jE _i) with allocations P(z _i=j)=w _j

and

$$w\sim{\cal D}\,(\delta _{{\rm 1}} ,\ldots,\,\delta _{k} )$$

follows a Dirichlet distribution. We also make the additional assumption that the δj’s are equal to 1, so that p(w) is a uniform distribution on the space described by (4). For the Poisson parameters λ _j, we take Gamma priors, so that

$$\lambda _{j} \sim{\rm Gamma}\,(a,\,b),\quad \quad j\,=\,{\rm 1},\ldots,k$$

with the ordering constraint

(5) $$\lambda _{{\rm 1}} \,\lt \,\lambda _{{\rm 2}} \,\lt \,\cdots \lt \,\lambda _{k} $$

to ensure that the components are identifiable. The ordering constraint is not necessary for the Monte Carlo algorithm to work. However, it does avoid the problem of label switching, as otherwise, any permutation of the indices {1, … ,k} will result in the same posterior distribution. The choice of objective priors is particularly difficult for finite mixture models, as common improper priors will lead to improper posteriors; see subsection 3.2.2 of Frühwirth-Schnatter (Reference Green2006) for a discussion. Note that a Dirichlet distribution as the prior on the weights is conjugate so we also have a Dirichlet posterior distribution for the weights but with updated parameters. In addition, we choose identical priors on the hyperparameters λ _j so the model is invariant to ordering.

There are k! ways to order k distinct objects. Because of the ordering constraint in equation (5), and the fact that the λs are i.i.d., the joint density of the collective λ is

$$p\,(\lambda \mid\alpha ,\,\beta ,\,k)\,=\,k\,!\,\,p\,(\lambda _{1} \mid\alpha ,\,\beta )\, \cdots \,p\,(\lambda _{k} \mid\alpha ,\,\beta )\,I_{{\lambda _{1} \,} {{\lt \,\lambda _{2} \,\lt \,\cdots\lt \lambda _{k} }}} \,(\lambda )$$

When k is fixed and known, the factorial term k! does not affect the MCMC algorithm as it can be absorbed into the normalising constant. However, in the variable k case, it must be noted, as it is a factor in the reversible jump acceptance probability. The joint density of all unknowns is

(6) $$\pi \,(w,\,\lambda ,\,z\mid D^{n} )\,\propto\,p\,(w\mid\delta )\,p\,(z\mid w)\,p\,(\lambda \mid\alpha ,\,\beta )\,L\,(D^{n} \mid\lambda ,\,z,\,E^{n} )$$

With the missing data formulation, the likelihood term L (D ⁿ ∣ λ, z, E ⁿ) can be written as

$$L\,(D^{n} \mid\lambda ,\,z,\,E)\,=\,\prod\limits_{i\,=\,1}^n {\left( {{{e^{{{\minus}\lambda _{{z_{i} }} E_{i} }} \,(\lambda _{{z_{i} }} E_{i} )^{{D_{i} }} } \over {D_{i} \,!\,}}} \right)} $$

and

$$p\,(z\mid w)\,=\,\prod\limits_{j\,=\,1}^k {w_{j}^{{n_{j} }} } $$

where n _j=#{i ∣ z _i=j} is the number of observations allocated to component j. The prior distributions are

$$\eqalignno{ & p\,(w\mid\delta )\,=\,{{\Gamma \,\left( {\Sigma _{{j\,=\,1}}^{k} \delta _{j} } \right)} \over {\prod_{{j\,=\,1}}^{k} \Gamma \,(\delta _{j} )}}\prod\limits_{j\,=\,1}^k {w_{j}^{{\delta _{j} \,{\minus}\,1}} } \cr & p\,(\lambda _{i} \mid\alpha ,\,\beta )\,=\,{{\beta ^{\alpha } } \over {\Gamma \,(\alpha )}}\lambda _{i}^{{\alpha \,{\minus}\,1}} e^{{{\minus}\beta \lambda _{i} }} $$

Using Bayes’ theorem, we have the following posterior conditional distributions:

$$\pi \,(\lambda _{j} )\,\propto\,\lambda _{j}^{{\alpha \,{\minus}\,1}} e^{{{\minus}\beta \lambda _{j} }} \,{\times}\,\lambda _{j}^{{\Sigma _{{i\,\mid\, z_{i} \,=\,j}} D_{i} }} e^{{{\minus}\lambda _{j} \Sigma _{{i\,\mid\, z_{i} \,=\,j}} E_{i} }} \,I_{{(\lambda _{{j\,{\minus}\,1,\,}} \lambda _{{j\,{\plus}\,1}} )}} \,(\lambda _{j} )$$

and

$$\pi \,(w\mid\delta ,\,z)\,\propto\,p(w\mid\delta )\,p\,(z\mid w)$$

so that

$$w\sim{\cal D}\,(\delta _{{\rm 1}} \,{\plus}\,n_{{\rm 1}} ,\ldots,\delta _{k} \,{\plus}\,n_{{_{k} }} )$$

where n _j=#{i|z _i=j}. For z, we update the allocations using

$$P\,(z_{i} \,=\,j)\,\propto\,w_{j} f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )\quad i\,=\,1,\ldots,n;\quad j\,=\,1,\ldots,k$$

so that

(7) $$p\,(z_{i} \,=\,j)\,=\,{{w_{j} f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )} \over {\Sigma _{{j\,=\,1}}^{k} \,w_{j} f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )}}$$

This follows from equation (2). The Gibbs algorithm for fixed k is then (Robert & Casella, Reference Robert and Casella1999)

Step 1: Simulate z _i from

$$p\,(z_{i} \,=\,j)\,\propto\,w_{j} f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )\ {\rm for}\ j\,=\,{\rm 1},\ldots,k$$

and compute $n_{j} ,\,n_{j} \bar{D}_{j} ,\,n_{j} \bar{E}_{j} $ from

$$n_{j} \,=\,\mathop{\sum}\limits_{i\mid z_{i} \,=\,j} {(1)} \quad n_{j} \bar{D}_{j} \,=\,\mathop{\sum}\limits_{i\mid z_{i} \,=\,j} {D_{i} } \quad n_{j} \bar{E}_{j} \,=\,\mathop{\sum}\limits_{i\mid z_{i} \,=\,j} {E_{i} } $$

Step 2: Simulate

$$\lambda _{j} \sim{\rm Gamma}\,(\alpha \,{\plus}\,n_{j} \bar{D}_{j} ,\,\beta \,{\plus}\,n_{j} \bar{E}_{j} )\,{\Bbb I}_{{(\lambda _{{j\,{\minus}\,1}} ,\,\lambda _{{j\,{\plus}\,{\rm 1}}} )}} \,(\lambda _{j} )\;{\rm for}\;j\,=\,{\rm 1},\ldots,k$$

Step 3: Simulate

$$w\sim{\cal D}\,(\delta \,{\plus}\,n_{{\rm 1}} ,\ldots,\delta \,{\plus}\,n_{k} )$$

3.1. RJMCMC

The reversible jump algorithm is an extension of the Metropolis–Hastings algorithm. We assume there is a countable collection of candidate models, indexed by $$M_{i} \in {\scr M}\,=\,\{ M_{1} ,\,M_{2} ,\ldots,M_{k} \} $$ . We further assume that for each model M _i, there exists an unknown parameter vector $\theta _{i} \in {\Bbb R}^{{n_{i} }} $ , where n ⁱ, the dimension of the parameter vector, can vary with i. Typically, we are interested in finding which models have the greatest posterior probabilities, and also to obtain estimates of the parameters. Thus, the unknowns in this modelling scenario will include the model index M _i, as well as the parameter vector θ _i. We assume that the models and corresponding parameter vectors have a joint density π(M _i, θ _i). The reversible jump algorithm constructs a reversible Markov chain on the state space $${\scr M}\,{\times}{\cup}_{{M_{i} \,\in \,{\scr M}}} \,{\Bbb R}^{{n_{i} }} $$ , which has π as its stationary distribution (Green, Reference Green and Richardson1995). In many instances, and in particular for Bayesian problems, this joint distribution is of the form

$$\pi \,(M_{i} ,\,\theta _{i} )\,=\,\pi \,(M_{i} ,\,\theta _{i} \mid X)\,\propto\,L(X\mid M_{i} ,\,\theta _{i} )\,p\,(M_{i} ,\,\theta _{i} )$$

where the prior on (M _i, θ _i) is often of the form

$$p\,(M_{i} ,\,\theta _{i} )\,=\,p(\theta _{i} \mid M_{i} )\,p\,(M_{i} )$$

with p(M _i) being the density of some counting distribution.

Suppose we are at model M _i and a move to model Mj is proposed, with probability r _ij. The corresponding move from θ _i to θ _j is achieved by using a deterministic transformation h _ij, such that

(8) $$(\theta _{j} ,\,v)\,=\,h_{{ij}} \,(\theta _{i} ,\,u)$$

where u and v are random variables introduced to ensure dimension matching necessary for reversibility. To ensure dimension matching, it is necessary that

$$\dim \,(\theta _{j} )\,{\plus}\,\dim \,(v)\,=\,\dim \,(\theta _{i} )\,{\plus}\,\dim \,(u)$$

For discussions about possible choices for the function h _ij, we refer the reader to Green (Reference Green and Richardson1995) and Brooks et al. (Reference Brooks, Giudici and Roberts2003b). Let

(9) $$A\,(\theta _{i} ,\,\theta _{j} )\,=\,{{\pi \,(M_{j} ,\,\theta _{j} )} \over {\pi \,(M_{i} ,\,\theta _{i} )}}\,{{q\,(v)} \over {q\,(u)}}{{r_{{ji}} } \over {r_{{ij}} }}\;\left| {{{\partial h_{{ij}} \,(\theta _{i} ,\,u)} \over {\partial \,(\theta _{i} ,\,u)}}} \right|$$

The acceptance probability for a proposed move from model (M _i, θ _i) to model (M _j, θ _j) is then

$${\rm min}\,\{ {\rm 1},\,A\,(\theta _{i} ,\,\theta _{j} )\} $$

where q(u) and q(v) are the respective proposal densities for u and v, and ∣ ∂h _ij (θ _i, u)/∂(θ _i, u) ∣ is the Jacobian of the transformation induced by h _ij. Green (Reference Green and Richardson1995) shows that the algorithm with acceptance probability given above simulates a Markov chain that is reversible and follows from the detailed balance equation

$$\pi \,(M_{i} ,\,\theta _{i} )\,q\,(u)\,r_{{ij}} \,=\,\pi (M_{j} ,\,\theta _{j} )\,q\,(v)\,r_{{ji}} \left| {{{\partial h_{{ij}} \,(\theta _{i} ,\,u)} \over {\partial \,(\theta _{i} ,\,u)}}} \right|$$

Detailed balance is necessary to ensure reversibility and is a sufficient condition for the existence of a unique stationary distribution. For the reverse move from model M _j to model M _i, it is easy to see that the transformation used is (θ _i, u)=h ⁻¹(θ _j, v) and the acceptance probability for such a move is

$$\min \,\left\{ {1,\,{{\pi \,(M_{i} ,\,\theta _{i} )} \over {\pi \,(M_{j} ,\,\theta _{j} }}\,{{q\,(u)} \over {q\,(v)}}\,{{r_{{ij}} } \over {r_{{ji}} }}\left| {{{\partial h_{{ij}} \,(\theta _{i} ,\,u)} \over {\partial \,(\theta _{i} ,\,u)}}} \right|^{{{\minus}1}} } \right\}\,=\,\min \,\{ 1,\,A\,(\theta _{i} ,\,\theta _{j} )^{{{\minus}1}} \} $$

For inference regarding which model has the greater posterior probability, we can base our analysis on a realisation of the Markov chain constructed above. The marginal posterior probability of model M _i is

$$\pi \,(M_{i} \mid X)\,=\,{{p\,(M_{i} )\,f\,(X\mid M_{i} )} \over {\Sigma _{{M_{j} \,\in \,{\scr M}}} \,p\,(M_{j} )\,f\,(X\mid M_{j} )}}$$

where

$$f\,(X\mid M_{i} )\,=\,\mathop{\int} {L\,(} X\mid M_{i} ,\,\theta _{i} )\,p\,(\theta _{i} \mid M_{i} )\,d\theta _{i} $$

is the marginal density of the data, which is obtained by integrating over the unknown parameters θ. In practice, we estimate π(M _i ∣ X) by counting the number of times the Markov chain visits model M _i in a single long run after reaching stationarity. These between-model moves described in this section are also augmented with within-model Gibbs updates as given in section 2.3 to update model parameters.

To assess convergence of the reversible jump algorithm, we use the method of Brooks & Giudici (Reference Brooks and Giudici1999) in which they propose to run I≥2 chains in parallel and base their convergence diagnostic on splitting the total variation not just between chains, but also between models. Their method was extended by Brooks et al. (Reference Brooks, Giudici and Philippe2003a) to include non-parametric techniques, including χ ² tests, Kolmogorov–Smirnov tests, and direct convergence rate estimation.

Brooks et al. (Reference Brooks, Giudici and Philippe2003a) suggest several methods for assessing convergence within the context of model selection problems. In particular, for reversible jump algorithms, we can have some idea of how fast the simulations approach stationarity by comparing the empirical stationary distribution with the observed model orders. They propose specific test statistics based on the χ ² distribution and also a Kolmogorov–Smirnov test for goodness of fit. The χ ² and Kolmogorov–Smirnov compare the stationary distribution of each chain and computes p-values for the computed test statistics. A critical value of 5% is used so that if the χ ² or Kolmogorov–Smirnov statistic is above this significance level there is no reason to reject the chains as not being from the same stationary distribution. See Brooks et al. (Reference Brooks, Giudici and Philippe2003a) for further details.

3.2. Reversible jump model selection

To update the model order, and thereby increase or decrease the number of components in the mixture, we use a combination of birth/death and split/merge moves as described below. We assume a uniform prior on the number of components k, so that

$$k\sim U\,\{ {\rm 1},\ldots,k_{{{\rm max}}} \} $$

where k _max is chosen to allow the algorithm to explore all feasible models. We set k _max=72, the number of groups, as under our hypothesis, this is the maximum number of components in the mixture. k=k _max only when the groups are all distinct. Setting k _max=72 will allow for direct comparison of the empirical Bayesian and the mixture model approach.

By introducing a prior on the number of components k, we extend the joint density (6) of all parameters. This yields

(10) $$\pi \,(k,\,w,\,z,\,\lambda \mid D^{n} )\,\propto\,p\,(k)\,p\,(w\mid\delta ,\,k)\,p\,(\lambda \mid\alpha ,\,\beta ,\,k)\,p\,(z\mid k)\,L\,(D^{n} \mid\lambda ,\,z)$$

Note that the densities of the other model parameters now depend on k. In sections 3.3 and 3.5, we describe in detail two algorithms that are used to simulate from this density. These algorithms are then combined with the fixed k updates of section 2.3 to simulate from the density in equation (10). Modelling mixtures with and without the Dirichlet process prior is considered by Green & Richardson (Reference Green and Richardson2001), who also consider the case of an unknown number of components. Alternatives to the reversible jump algorithm exist in this context. For example, Dellaportas & Karlis (Reference Dellaportas and Karlis2001) develop a semi-parametric sample-based method to approximate a mixing density g(θ), based on the method of moments.

3.3. Split and merge moves

Note that the joint density in equation (10) now depends on k. We use the split/merge method of Dellaportas et al. (Reference Dellaportas, Forster and Ntzoufras1997) and Viallefont et al. (Reference Viallefont, Richardson and Green2002). Suppose we are at a configuration with k components, with

$$\theta _{k} \,=\,\{ (\lambda _{{\rm 1}} ,\,w_{{\rm 1}} ),\ldots,\,(\lambda _{k} ,\,w_{k} )\} $$

and suppose a move to increase the number of components is proposed. We select uniformly one of the current k components to be split. Suppose the jth component (λ _j, w _j) is selected to be split into two components $$(\lambda _{{j_{1} }} ,\,w_{{j_{1} }} )$$ and $$(\lambda _{{j_{2} }} ,\,w_{{j_{2} }} )$$ such that j ₁=j and j ₂=j+1, the components originally numbered j+1, … , k are then renumbered j+2, … , k+1. The split is also designed so that the first two moments of the split component remain the same as the original component. Thus, we simulate u1 and u2 from densities defined on the interval [0, 1]. Usually, we use β densities and set

$$\eqalignno{ & w_{{j_{{\rm 1}} }} \,=\,w_{j} u_{{\rm 1}} \cr & w_{{j_{{\rm 2}} }} \,=\,w_{j} \,({\rm 1}\,{\minus}\,u_{{\rm 1}} ) \cr & \lambda _{{j_{{\rm 1}} }} \,=\,\lambda _{j} u_{{\rm 2}} \cr & \lambda _{{j_{{\rm 2}} }} \,=\,\lambda _{j} \,({\rm 1}\,{\minus}\,u_{{\rm 1}} u_{{\rm 2}} )/({\rm 1}\,{\minus}\,u_{{\rm 1}} ) $$

Other choices for splitting and merging components are described in Viallefont et al. (Reference Viallefont, Richardson and Green2002). The proposed parameter is then

$$\eqalignno{ \theta _{{k{\plus}{\rm 1}}} = & \{ (\lambda _{{\rm 1}} ,\,w_{{\rm 1}} ),\ldots,\,(\lambda _{j} _{{{\minus}{\rm 1}}} ,\,w_{{j{\minus}{\rm 1}}} ),\,(\lambda _{{j_{{\rm 1}} }} ,\,w_{{j_{{\rm 1}} }} ),\,(\lambda _{{j_{{\rm 2}} }} ,\,w_{{j_{{\rm 2}} }} ) \cr & (\lambda _{{j{\plus}{\rm 1}}} ,\,w_{{j{\plus}{\rm 1}}} ),\ldots,\,(\lambda _{k} ,\,w_{k} )\} $$

If the ordering constraint in equation (5) is not satisfied, then the move is rejected immediately, as the reverse move in which we merge two adjacent components would not be possible. We can compute the Jacobian for this transformation as

(11) $$\left| {{{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{k} ,\,u_{1} ,\,u_{2} )}}} \right|=\left| {{{\partial \,(w_{{j_{1} }} ,\,w_{{j_{2} }} ,\,\lambda _{{j_{1} }} ,\,\lambda _{{j_{2} }} )} \over {\partial \,(w_{j} ,\,\lambda _{j} ,\,u_{1} ,\,u_{2} )}}} \right|={{\lambda _{j} w_{j} } \over {1{\minus}u_{1} }}$$

For the reverse move, we select a pair of adjacent components j ₁ and j ₂. Combining them to get a new component labelled j, we set

$$w_{j} =w_{{j_{1} }} {\plus}w_{{j_{2} }} ,\quad \lambda _{j} ={{w_{{j_{1} }} \lambda _{{j_{1} }} {\plus}w_{{j_{2} }} \lambda _{{j_{2} }} } \over {w_{{j_{1} }} {\plus}w_{{j_{2} }} }}$$

by keeping the first two moments of the proposed and current configuration constant. We then sample a new set of allocation variables according to equation (7). We also keep track of the probability of each allocation, so that p _a(z) represents the probability of a given allocation. To compute p _a(z), we first simulate z _i using equation (7). For each i, the probability of that allocation is given by

$$p_{a} \,(z_{i} )={{w_{{z_{i} }} f\,(D_{i} \mid\lambda _{{z_{i} }} ,\,E_{i} )} \over {\Sigma _{{j=1}}^{k} w_{j} f\,(D_{i} \mid\lambda _{j} ,\,E_{i} )}}$$

Finally, we compute the probability of all allocations by

$$p_{a} \,(z)=\prod\limits_{i=1}^n {p_{a} \,(z_{i} )} $$

3.4. Acceptance probability

The acceptance probability of a move of type (k, θ _k) ⇒ (k′, θ _k′) is then min{1, A _{k, k′}}, where

$$\eqalignno{ A_{{k,\,k'}} =\, & {{\pi \,(k',\,\theta _{{k'}} )} \over {\pi \,(k,\,\theta _{k} )}}{\times}{{p\,(k'\Rightarrow k)} \over {p\,(k\Rightarrow k')}}{\times}{1 \over {q\,(u_{1} )q\,(u_{2} )}}{\times}\left| {{{\partial \theta _{{k'}} } \over {\partial \,(\theta _{k} ,\,u_{1} ,\,u_{2} )}}} \right| \cr A_{{k,\,k'}} = & \,{{p\,(k')p\,(w'\mid\delta ,\,k')p\,(\lambda '\mid\alpha ,\,\beta ,\,k')L\,(D^{n} \mid\lambda ',\,z')} \over {p\,(k)p\,(w\mid\delta ,\,k)p\,(\lambda \mid\alpha ,\,\beta ,\,k)L\,(D^{n} \mid\lambda ,\,z)}}\cr &{\times} {{p(z'\mid w',\,k{\plus}1)/p_{a} \,(z')} \over {p\,(z\mid w,\,k)/p_{a} \,(z)}}{\times}{{p\,(k'\Rightarrow k)} \over {p\,(k\Rightarrow k')}}{1 \over {q\,(u_{1} )q\,(u_{2} )}}\left| {{{\partial \theta _{{k'}} } \over {\partial \,(\theta _{k} ,\,u_{1} ,\,u_{2} )}}} \right| $$

with k′=k+1 this becomes

$$\eqalignno{ A_{{k,\,k{\plus}1}} = & \,{{p\,(k{\plus}1)} \over {p\,(k)}}{\times}{{p\,(w'\mid\delta ,\,k{\plus}1)} \over {p\,(w\mid\delta ,\,k)}}{\times}{{p\,(\lambda '\mid\alpha ,\,\beta ,\,k{\plus}1)} \over {p\,(\lambda \mid\alpha ,\,\beta ,\,k)}}\cr &{\times} {{L\,(D^{n} \mid\lambda ',\,w')} \over {L\,(D^{n} \mid\lambda ,\,w)}}{\times}{{p\,(k{\plus}1\Rightarrow k)} \over {p\,(k\Rightarrow k{\plus}1)}}{\times}{1 \over {q\,(u_{1} )q\,(u_{2} )}}\left| {{{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{{k}},\,{\rm u}_{1} ,\,{\rm u}_{2} )}}} \right| $$

Now with a uniform prior on the number of components k and the weights w

$$\eqalignno{ A_{{k,\,k{\plus}1}} = & {{\Gamma \,(k{\plus}1)} \over {\Gamma \,(k)}}{\times}{{(k{\plus}1)p\,(\lambda _{{j_{1} }} \mid\alpha ,\,\beta )p\,(\lambda _{{j_{2} }} \mid\alpha ,\,\beta )} \over {p\,(\lambda _{j} \mid\alpha ,\,\beta )}} \cr & {\times} {{p\,(z'\mid w',\,k{\plus}1)/p_{a} \,(z')} \over {p\,(z\mid w,\,k)/p_{a} \,(z)}}{\times}{{L\,(D^{n} \mid\lambda ',\,z')} \over {L\,(D^{n} \mid\lambda ,\,z)}}{\times}{{m_{{k{\plus}1}} } \over {s_{k} }} \cr & {\times}{1 \over {q\,(u_{1} )q\,(u_{2} )}}{\times}\left| {{{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{{k}},\,{\rm u}_{1} ,\,{\rm u}_{2} )}}} \right| $$

where the ratio of Gamma terms comes from the ratio of the prior distributions on $w\hskip0.6pt'$ and w and

$${{p\,(k{\plus}1\Rightarrow k)} \over {p\,(k\Rightarrow k{\plus}1)}}={{m_{{k{\plus}1}} /(k{\plus}1{\minus}1)} \over {s_{k} /k}}={{m_{{k{\plus}1}} } \over {s_{k} }}$$

3.5. Birth and death moves

Suppose we are now at model M _k with k components, say

(12) $$\theta _{k} =\left\{ {(\lambda _{{\rm 1}} ,\,w_{{\rm 1}} ),\ldots,(\lambda _{k} ,\,w_{k} )} \right\}$$

If a move is proposed to increase the number of components by one, then we simulate

$$\tilde{w}\sim{\rm Beta}({\rm 1},\,k)\ {\rm and}\ \tilde{\lambda }\sim{\rm Gamma}(a,\,b)$$

independently. The proposed new component will then have weight $\tilde{w}$ and the other weights are then scaled by a factor of $({\rm 1}{\minus}\tilde{w})$ , so that the sum of the weights remain equal to 1. The corresponding Poisson parameter for the proposed component in $\tilde{\lambda }$ . Note that $\tilde{\lambda }$ is sampled from its prior distribution. The proposed component is then

(13) $$\theta _{{k{\plus}{\rm 1}}} =(\lambda _{{\rm 1}} ,\,w_{{\rm 1}} /({\rm 1}{\minus}\tilde{w})),\ldots,(\lambda _{{k,\,w_{k} }} /({\rm 1}{\minus}\tilde{w})),\,(\tilde{\lambda },\,\tilde{w})\} $$

Using this proposed value of θ _k+1, we also simulate proposed values for the allocations $z\hskip0.1pt'$ with model k+1. Using the general form of the reversible jump acceptance probability, see, for example, Green (Reference Green and Richardson1995), the probability of changing the number of components to k+1 is then min{1, A _{k, k+1}}, where

$$A_{{k,\,k{\plus}1}} ={{\pi \,(k{\plus}1,\,\theta _{{k{\plus}1}} )} \over {\pi (k,\,\theta _{k} )}}{\times}{{p\,(k{\plus}1\Rightarrow k)} \over {p\,(k\Rightarrow k{\plus}1)}}{\times}{1 \over {q\,(\tilde{w})q\,(\tilde{\lambda })}}{\times}\left| {{{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{k} ,\,\tilde{w},\,\tilde{\lambda })}}} \right|$$

Making the necessary substitutions yield

(14) $$\eqalignno{ A_{{k,\,k{\plus}1}} = \,& {{p\,(k{\plus}1)p\,(w'\mid\delta ,\,k{\plus}1)p\,(\lambda '\mid\alpha ,\,\beta ,\,k{\plus}1)L\,(D^{n} \mid\lambda ',\,z')} \over {p\,(k)p\,(w\mid\delta ,\,k)p\,(\lambda \mid\alpha ,\,\beta ,\,k)L\,(D^{n} \mid\lambda ,\,z)}}\cr & {\times} {{p\,(z'\mid w',\,k{\plus}1)/p_{a} \,(z')} \over {p\,(z\mid w,\,k{\plus}1)/p_{a} \,(z)}}{\times}{{p\,(k{\plus}1\Rightarrow k)} \over {p\,(k\Rightarrow k{\plus}1)}}{\times}{1 \over {q\,(\tilde{w})q\,(\tilde{\lambda })}} \cr & {\times}\left| {{{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{k} ,\,\tilde{w},\,\tilde{\lambda })}}} \right| $$

Using equations (12) and (13) we then have the Jacobian

$${{\partial \theta _{{k{\plus}1}} } \over {\partial \,(\theta _{k} ,\,\tilde{w},\,\tilde{\lambda })}}=(1{\minus}\tilde{w})^{{k{\minus}1}} $$

If we denote the probability of a birth when there are k components by b _k, and the probability of a death by d _k, with b _k+d _k=1, then

$${{p\,(k{\plus}1\Rightarrow k)} \over {p\,(k\Rightarrow k{\plus}1)}}={{d_{{k{\plus}1}} /(k{\plus}1)} \over {b_{k} }}$$

as for the move to be reversible we would then be able to kill k+1 components in the new model, each with equal probability. Substituting these values in equation (14), the ratio A _k,k+1 reduces to

$$\eqalignno{ A_{{k,\,k{\plus}1}} = & {{\Gamma \,(k{\plus}1)} \over {\Gamma \,(k)}}{\times}(k{\plus}1)p\,(\tilde{\lambda }){\times}{{L\,(D^{n} \mid\lambda ',\,z')} \over {L\,(D^{n} \mid\lambda ,\,z)}}{\times}{{p\,(z'\mid w',\,k{\plus}1)/p_{a} \,(z')} \over {p\,(z\mid w,\,k)/p_{a} \,(z)}} \cr {\times} & {{d_{{k{\plus}1}} /(k{\plus}1)} \over {b_{k} }}{1 \over {q\,(\tilde{w})q\,(\tilde{\lambda })}}{\times}(1{\minus}\tilde{w})^{{k{\minus}1}} $$

which on substituting $q(\tilde{\lambda })=p(\tilde{\lambda })$ and $q\,(\tilde{w})=k\,({\rm 1}{\minus}\tilde{w})^{{k{\minus}{\rm 1}}} $ further reduces to

(15) $$A_{{k,\,k{\plus}1}} ={{p\,(z'\mid w',\,k{\plus}1)/p_{a} \,(z')} \over {p\,(z\mid w,\,k)/p_{a} \,(z)}}{{L\,(D^{n} \mid\lambda ',\,z')} \over {L\,(D^{n} \mid\lambda ,\,z)}}{\times}{{d_{{k{\plus}1}} } \over {b_{k} }}$$

For a proposed death move, the acceptance probability is then

$${\rm min}\{ 1,\,A_{{k,\,k{\plus}1}}^{{{\minus}1}} \} $$

Although the algorithm simulates new values of the allocations when proposing to move, it is not necessary to carry the allocations along. For between-model moves, we could replace the missing data formulation by noting that

$${{p\,(z\mid w,\,k)} \over {p_{a} \,(z)}}L\,(D^{n} \mid\lambda ,\,z)=L\,(D^{n} \mid\lambda ,\,w)$$