
\documentclass{article}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{amssymb}
\usepackage{amsmath}

\setcounter{MaxMatrixCols}{10}
%TCIDATA{OutputFilter=LATEX.DLL}
%TCIDATA{Version=5.00.0.2606}
%TCIDATA{<META NAME="SaveForMode" CONTENT="1">}
%TCIDATA{BibliographyScheme=Manual}
%TCIDATA{Created=Wednesday, December 28, 2005 18:33:59}
%TCIDATA{LastRevised=Sunday, January 01, 2006 21:34:32}
%TCIDATA{<META NAME="GraphicsSave" CONTENT="32">}
%TCIDATA{<META NAME="DocumentShell" CONTENT="Standard LaTeX\Blank - Standard LaTeX Article">}
%TCIDATA{CSTFile=40 LaTeX article.cst}

\newtheorem{theorem}{Theorem}
\newtheorem{acknowledgement}[theorem]{Acknowledgement}
\newtheorem{algorithm}[theorem]{Algorithm}
\newtheorem{axiom}[theorem]{Axiom}
\newtheorem{case}[theorem]{Case}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{conclusion}[theorem]{Conclusion}
\newtheorem{condition}[theorem]{Condition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{criterion}[theorem]{Criterion}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{exercise}[theorem]{Exercise}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{notation}[theorem]{Notation}
\newtheorem{problem}[theorem]{Problem}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{solution}[theorem]{Solution}
\newtheorem{summary}[theorem]{Summary}
\newenvironment{proof}[1][Proof]{\noindent\textbf{#1.} }{\ \rule{0.5em}{0.5em}}
\input{tcilatex}

\begin{document}


2.11 Monte Carlo Estimates of Model Probabilities\newline
\newline
Here a Monte Carlo method for estimating posterior model probabilities based
on independent MCMC sampling of two (or more) different models; this is a
modified version of the approach suggested in Congdon (2006). Let $\theta
=(\theta _{1},...,\theta _{J}$) denote the parameter set over all J models,
with dimension (d$_{1}$,..d$_{J}$). Assume a model indicator m $\in $%
(1,\dots ,J), such that given m=j, $\theta _{j}$ defines the likelihood for
Y=(y$_{1}$,...y$_{N}$) and Y is independent of parameters of other models $%
\theta _{k},k\neq j$ (Godsill, 2001; Carlin \&\ Chib, 1995). This means that
the marginal likelihood given m=j is $\qquad \qquad P(Y|m=j)=\int P(Y|\theta
,m=j)P(\theta |m=j)d\theta =\int P(Y|\theta _{j},m=j)P(\theta |m=j)d\theta
_{j}$\newline
\linebreak This framework replicates that of section 2.1 in Carlin \&\ Chib
(1995). As these authors mention, the form of the cross-model priors P($%
\theta _{k}|m=j,j\neq k)$ within the product P($\theta |m=j)=\underset{k}{%
\prod }$P($\theta _{k}|m=j)$ is arbitrary, though proper densities are
required in order that P($\theta |m=j)$ integrates to 1. Here the assumption
is made that P($\theta _{k}|m=j_{1})=$P($\theta _{k}|m=j_{2})$ for all k$%
\neq j_{1},$ k $\neq j_{2},$ so that there will be J cross-model priors. So
in a three model situation P($\theta _{3}|m=1)=P(\theta _{3}|m=2)=g_{3},$P($%
\theta _{2}|m=1)=P(\theta _{2}|m=3)=g_{2},$ and P($\theta _{1}|m=2)=P(\theta
_{1}|m=3)=g_{1}.$ This choice of cross-model prior contrasts with the flat
option for P($\theta _{k}|m=j,k\neq j)$ considered by Congdon (2006).\newline
\linebreak With this framework, one may use the output \{$\theta
_{j}^{(t)},t=1,..,T;j=1,..,J\}$ from the J models to estimate $P$(m=j{%
\TEXTsymbol{\vert}}Y). Note that this estimate needs a sample of the same
length (say T iterations) from the posteriors P($\theta _{j}|Y,m=j)$ of all
J models under consideration. Such samples might be obtained by running
models in parallel or by running them separately and then pooling the output
(e.g. in a spreadsheet). This is conceptually distinct from product space
search algorithms, such as that of Carlin and Chib (1995), when the model j
parameters are only updated when model j is visited. To see how a Monte
Carlo posterior model probability is estimated from such output, first write

$\qquad P$(m=j{\TEXTsymbol{\vert}}Y)=$\int P$(m=j,$\theta $\TEXTsymbol{\vert}%
Y$)d\theta =\int P$(m=j{\TEXTsymbol{\vert}}Y,$\theta $)P($\theta ${%
\TEXTsymbol{\vert}}Y)d$\theta $. \newline
Hence a Monte Carlo estimate of P(m=j{\TEXTsymbol{\vert}}Y) is obtainable as

$\qquad w_{j}$ = $\overset{T}{\underset{t=1}{\sum }}P$(m=j{\TEXTsymbol{\vert}%
}Y,$\theta $$^{(t)}$)/T \newline
where \{$\theta $$^{(t)}$= ($\theta _{1}^{(t)}$,\dots ,$\theta _{2}^{(t)}$%
,\dots ,$\theta _{J}^{(t)}$), t=1,T\} are T samples of parameters in all
models. For a particular iteration let

$w_{j}^{(t)}=P(m=j{|}Y,\theta ^{(t)})=\frac{P(m=j,Y,\theta ^{(t)})}{%
P(Y,\theta ^{(t)})}=\frac{P(Y{|}m=j,\theta ^{(t)})P(\theta ^{(t)}{|}%
m=j)P(m=j)}{P(Y,\theta ^{(t)})}\qquad (2.15)$\newline
The numerator in (2.15) contains the term $P(\theta |m=j)$=$P(\theta
_{1},\theta _{2},..\theta _{j,}...\theta _{J}|m=j).$ From above the
cross-model prior is arbitrary and the simplifying assumption

$\ P(\theta _{h}{|}m=k)=g_{h}$ \ \ \ \ \ \ ($all$ $h\neq k)$\newline
is made, where g$_{h}$ is a proper density. So \newline
$P(\theta _{h}{|}m=1)=P(\theta _{h}{|}m=2)=..=P(\theta _{h}{|}%
m=h-1)=P(\theta _{h}{|}m=h+1)..P(\theta _{h}{|}m=J)=g_{h}$. \newline
and there are J cross model priors \{g$_{1}$,..,g$_{J}$\}. As in the model
choice procedures of Gelfand and Dey (1994) and Carlin \&\ Chib (1995), one
might estimate g$_{h}$ as an approximation to P($\theta _{h}|m=h,Y),$ namely
the posterior density of $\theta _{h}$ given Y and m=h. It follows that

$P(\theta |m=j)=P(\theta _{j}{|}m=j)\underset{h\neq j}{\overset{J}{\prod }}%
P(\theta _{h}{|}m=j)=P(\theta _{j}{|}%
m=j)[g_{1}g_{2}..g_{j-1}g_{j+1}..g_{J}]. $\linebreak Then

$\qquad w_{j}^{(t)}=\frac{P(Y{|}m=j,\theta _{j}^{(t)})P(\theta _{j}^{(t)}{|}%
m=j)[\underset{h\neq j}{\Pi }g\overset{(t)}{_{h}}]P(m=j)}{P(Y,\theta ^{(t)})}%
\qquad \qquad \qquad (2.16)$\newline
The denominator in (2.16) can be written

\qquad P(Y, $\theta $$^{(t)}$)=$\overset{J}{\underset{k=1}{\sum }}$P(Y, $%
\theta $$^{(t)}$, m=k)

\qquad \qquad \qquad =$\overset{J}{\underset{k=1}{\sum }}$\{P(Y{\TEXTsymbol{%
\vert}} $\theta $$^{(t)},m=k$) P($\theta _{k}^{(t)}{|}m=k$)$[\underset{h\neq
k}{\Pi }g\overset{(t)}{_{h}}]$P(m=k)\}

\qquad \qquad \qquad = $\overset{J}{\underset{k=1}{\sum }}$P(Y{\TEXTsymbol{%
\vert}}$\theta _{k}^{(t)},m=k$) P($\theta _{k}^{(t)}{|}m=k$)$[\underset{%
h\neq k}{\Pi }g\overset{(t)}{_{h}}]$P(m=k).\newline
Then

$\qquad w_{j}^{(t)}=\frac{P(Y{|}m=j,\theta _{j}^{(t)})P(\theta _{j}^{(t)}{|}%
m=j)[\underset{h\neq j}{\Pi }g\overset{(t)}{_{h}}]P(m=j)}{\overset{J}{%
\underset{k=1}{\sum }}P(Y{|}\theta _{k}^{(t)},m=k)P(\theta _{k}^{(t)}{|}m=k)[%
\underset{h\neq k}{\Pi }g\overset{(t)}{_{h}}]P(m=k)}$\newline
One may divide through by the product of the J cross-model priors [g$_{1}$,g$%
_{2},$...g$_{J}$] giving

$\qquad w_{j}^{(t)}=\frac{[\frac{P(Y{|}m=j,\theta _{j}^{(t)})P(\theta
_{j}^{(t)}{|}m=j)P(m=j)}{g\overset{(t)}{_{j}}}]}{\overset{J}{\underset{k=1}{%
\sum [}}\frac{P(Y{|}\theta _{k}^{(t)},m=k)P(\theta _{k}^{(t)}{|}m=k)P(m=k)}{g%
\overset{(t)}{_{k}}}]}$\newline
For example when J=2, one obtains

$w_{1}^{(t)}=\frac{[\frac{P(Y{|}m=1,\theta _{1}^{(t)})P(\theta _{1}^{(t)}{|}%
m=1)P(m=1)}{g\overset{(t)}{_{1}}}]}{[\frac{P(Y{|}m=1,\theta
_{1}^{(t)})P(\theta _{1}^{(t)}{|}m=1)P(m=1)}{g\overset{(t)}{_{1}}}]+[\frac{%
P(Y{|}m=2,\theta _{2}^{(t)})P(\theta _{2}^{(t)}{|}m=2)P(m=2)}{g\overset{(t)}{%
_{2}}}]}$\linebreak \linebreak \newline
While the above development builds on the framework of Carlin \&\ Chib
(1995) there are two main differences. First, the cross-model priors do not
function here as linking densities in a model search algorithm since there
is no switching between models. Instead all models are sampled from at each
iteration. This avoids the problems involved in tuning prior model
probabilities or 'jump' proposal densities in product space algorithms to
ensure that models are visited sufficiently often (Green \&\ O'Hagan, 1998;
Friel \&\ Pettitt, 2006). The second main difference follows from the first.
Switching between models under product space search implies a binary form of
model averaging: if model 1 is selected at several successive iterations
then the averaged parameters at these iterations are in fact the parameter
values sampled from one model only. Under the approach here averaging is on
the basis of continuous quantities, namely the $w_{k}^{(t)}$obtained for all
models, so some weight in the average at each iteration is given to inferior
models. This may be important when models are closely competing. The weights
are obtained by considering the relativities in log[P(Y,$\theta ,m=k)],$
namely the log of the product of the likelihood, the own model prior P($%
\theta _{k}|m=k)$ and the product of cross model priors . Thus at iteration
t, one obtains $q_{k}^{(t)}=\log \{P(Y{|}\theta _{k}^{(t)},m=k)P(\theta
_{k}^{(t)}{|}m=k)[\underset{a\neq k}{\Pi }g_{a}]P(m=k)\}$, and deviations$%
\Delta q_{k}^{(t)}=q_{k}^{(t)}-\underset{k}{max}(q_{k}^{(t)})$, with w$%
_{k}^{(t)}$ obtained by exponentiating:

\qquad \qquad $w_{k}^{(t)}$=$\frac{exp(\Delta q_{k}^{(t)})}{\underset{k}{%
\sum }exp(\Delta q_{k}^{(t)})}$ \newline
\newline
An indication of the uncertainty attaching to the estimated model weights
can be obtained by batch sampling (Lewis and Raftery, 1997). This involves
dividing the pooled output over T iterations from the J models into B
batches of equal size. A \ 95\% interval for the formal weight
approximations $\overset{-}{w}_{j}$ is based on an approximate t density
with B-1 degrees of freedom for the batch means\ $\overset{-}{w}%
_{jb},b=1,...,B$.\newline
\newline
This approach is illustrated with the radiata pine data considered earlier.
Estimates of the cross-model priors g$_{j,m}$ for parameters m and models j
are based on a preliminary run. A lognormal approximation is used for the
density of the precisions $\tau _{j}$ in equation (2.10), and univariate
normal approximations for the densities of $\alpha _{j},$ and $\beta _{j}.$
The relevant parameter values (posterior means and standard deviations from
the preliminary run) are in g.mn[] and g.se[] in the code below. With prior
probabilities Pr(m=1)=0.9995 and Pr(m=2)=0.0005, we find the posterior
probability of model 1 is estimated as 0.2953, similar to the value of
0.2914 reported by Green and O'Hagan (1998, p 8). The code is\newline
model \{for (i in 1:42) \{ x[i] \TEXTsymbol{<}- X[i]-mean(X[]); w[i] 
\TEXTsymbol{<}- W[i]-mean(W[])

for (j in 1:2) \{Y[i,j] \TEXTsymbol{<}- y[i]; Y[i,j] \symbol{126}%
dnorm(mu[i,j],tau[j])

\# log LKD contribution

LL[i,j] \TEXTsymbol{<}- -0.5*tau[j]*pow(y[i]-mu[i,j],2) +
0.5*log(tau[j]/6.2832) \}

\# means for M1 and M2

mu[i,1] \TEXTsymbol{<}- a[1]+ b[1]*x[i]; mu[i,2] \TEXTsymbol{<}- a[2] +
b[2]*w[i]\}

\# scaled likelihoods at each iteration

for (k in 1:2) \{ L[k] \TEXTsymbol{<}-
sum(LL[,k])+sum(logPr[1:3,k])+log(PrMod[k])-sum(g[1:3,k])

SL[k] \TEXTsymbol{<}- max(L[k]-maxL,-500); expSL[k] \TEXTsymbol{<}-
exp(SL[k]); w[k] \TEXTsymbol{<}- expSL[k]/sum(expSL[])\}

\# maximum of scaled likelihoods

maxL \TEXTsymbol{<}- ranked(L[],2)

for (j in 1:3) \{ for (k in 1:2) \{gT[j,k] \TEXTsymbol{<}-
1/pow(g.se[j,k],2)\}\}

\# priors

for (j in 1:2) \{

a[j] \symbol{126}dnorm(M[1],T[1]); b[j] \symbol{126}dnorm(M[2],T[2]) ;
tau[j] \symbol{126}dgamma(3,180000); logtau[j] \TEXTsymbol{<}- log(tau[j]);

\# probs on sampled parameters in terms of own model prior and cross-model
priors

logPr[1,j] \TEXTsymbol{<}- 0.5*log(T[1]/6.28)-0.5*T[1]*pow(a[j]-M[1],2)

logPr[2,j] \TEXTsymbol{<}- 0.5*log(T[2]/6.28)-0.5*T[2]*pow(b[j]-M[2],2)

logPr[3,j] \TEXTsymbol{<}-
log(pow(tau[j],r-1)*exp(-m*tau[j])*pow(m,r))-loggam(r)

g[1,j] \TEXTsymbol{<}-
0.5*log(gT[1,j]/6.28)-0.5*gT[1,j]*pow(a[j]-g.mn[1,j],2)

g[2,j] \TEXTsymbol{<}-
0.5*log(gT[2,j]/6.28)-0.5*gT[2,j]*pow(b[j]-g.mn[2,j],2)

g[3,j] \TEXTsymbol{<}-
0.5*log(gT[3,j]/6.28)-0.5*gT[3,j]*pow(log(tau[j])-g.mn[3,j],2)-log(tau[j])\}%
\}\newline
with initial value files list(a=c(3000,3000),b=c(200,200),tau=c(0.001,0.001))%
\newline
and list(a=c(3000,3000),b=c(100,100),tau=c(1,1)). Excluding observation
listings the data file is\
list(r=3,m=180000,T=c(0.000001,0.0001),M=c(3000,185),g.mn=structure(.Data=c(2991.6,2991.5,184.67,183.07,-11.612,-11.242),%
\newline
.Dim=c(3,2)),g.se=structure(.Data=c(51.72,43.148,11.46,9.4197,0.21223,0.21276),.Dim=c(3,2)),%
\newline
PrMod=c(0.9995,0.0005),y=c(..),X=c(..),W=c(..))\newline
\newline
For models with random effects b$_{i}$, let priors on fixed effects be
denoted $\eta $ and hyperparameters for the random effects be denoted $\Phi .
$ Then with $\theta =(\eta ,\Phi ),$ the posterior model probability is
estimated as $P$(m=k{\TEXTsymbol{\vert}}Y)= $\overset{-}{w}_{k}$= $\Sigma
w_{j}/T$ , where

\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $w_{j}^{(t)}=\frac{P(Y{|}m=j,\theta
_{j}^{(t)})P(\theta _{j}^{(t)}{|}m=j)[\underset{h\neq j}{\Pi }g\overset{(t)}{%
_{h}}]P(m=j)}{P(Y,\theta ^{(t)})}$\newline
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ P(Y\TEXTsymbol{\vert}m=k,$\theta )=$ $%
\int \overset{}{\underset{}{\overset{N}{\underset{i=1}{\Pi }}P(\text{y}_{%
\text{i}}|\beta ,b_{i})}}P(b_{i}|\Phi )db_{i}$ \newline
and the cross-model priors g$_{h}$ =P($\theta _{j}|m=h)$, h$\neq j,$ involve
(proper density) approximations for P$(\theta _{j}|Y,m=j)=$P$(\eta _{j},\Phi
_{j}|Y,m=j).$ In conjugate models (e.g. Poisson models with gamma random
effects) the integrated likelihoods P(Y\TEXTsymbol{\vert}m=k,$\theta )$ are
available analytically, whereas in general linear mixed models devices such
as numerical integration (e.g. by Simpson's rule) (Aitkin and Alf\`{o},
2003) or importance sampling (Geweke, 1989) are required, and are applied at
each iteration to obtain P(Y\TEXTsymbol{\vert}m=k,$\theta ^{(t)})$. For
example consider a scalar random effect in a Poisson lognormal model with y$%
_{i}$ \symbol{126}Po(E$_{i}\mu _{i})$ where E$_{i}$ are expected events, and
log($\mu _{i})$=X$_{i}\beta +b_{i},$ where b$_{i}\backsim $ N(0,1/$\tau
_{b}).$ Then\ an appropriate range for Simpson's integration over b$_{i}$
may be obtained from an earlier MCMC\ run.\newline

\end{document}
