In a previous post I looked at how to calculate marginal likelihoods for a simple model in R. For those of you interested in Julia, here is the example, ported to Julia. I haven't included the model switch integration - I leave that as an exercise, but feel free to post a comment if you want it posted here. The following, in the form of an IJulia notebook, can be downloaded as a gist.
To recap, under the Bayesian paradigm, inference is based on the posterior probability over the parameters of interest. It's helpful to think of our inferences being conditional on a given model, $M$ with a parameter vector $\theta \in \Theta$. Given a dataset, $D$, and a model, the posterior distribution of the parameter values is given by Bayes' theorem.
\begin{align} Pr(\theta|D,M) = \frac{Pr(D|\theta,M)Pr(\theta|M)}{Pr(D|M)} \end{align}
$Pr(D|\theta,M)$ is the likelihood function, $Pr(\theta|M)$ is the prior probability , and $Pr(D|M)$ is known as the marginal likelihood, predictive probability, or evidence. $Pr(D|M)$ is a normalising constant that ensures that $Pr(\theta|D,M)$ is a probability.
\begin{align} Pr(D|M) = \int_{\Theta}\Pr(D|\theta,M)Pr(\theta|M) d\theta \end{align}
However, one often wants to compare the fit of different models. As a function of the model $M$, the marginal likelihood can be interpreted sa the likelihood of the model $M$ given the data $D$. Hence, to choose between several models, one simply chooses the one with the highest marginal likelihood. When comparing two models,say $M_0$ and $M_1$, a ratio of marginal likelihoods, known as the Bayes factor, $K$, is usually defined:
\begin{align} K_{01} = \frac{Pr(D|M_1)}{Pr(D|M_0)} \end{align}
Interpretations of the Bayes factor have been provided by Jeffreys and Kass and Raftery.
Unfortunately, in most cases, it is not possible to obtain an exact solution of the marginal likelihood. A number of approaches have been described to obtain an approximate numerical estimate of the marginal likelihood; here, I illustrate two approaches based on tempering.
A simple example
Before I describe what tempering means in this context, let's consider a simple example, for which there is an analytical solution for the marginal likelihood. Consider the problem of fitting a set of $n=100$ exponential random variables, $X$ with parameter $\lambda=3$.
We can generate these in Julia as follows.
using Distributions
srand(1)
lambda=3.0
n=100
x=rand(Exponential(1.0/lambda),n);
However, as I want to compare with my R code, I'll use the same random data as I did in my R code.
x=[0.251727277709448, 0.393880926369035, 0.0485689089012643, 0.0465984206228326,
0.145356208593058, 0.964989512488025, 0.409854017804964, 0.179894279999038,
0.318855831232818, 0.0490153301679842, 0.463578376270143, 0.254009951823358,
0.41253451691161, 1.47464473922761, 0.351514389103556, 0.34508131534358,
0.625345057469416, 0.218248879071325, 0.11231115879491, 0.196159907151014,
0.788171751007882, 0.213964196077238, 0.0980401292132835, 0.18862184152628,
0.0353575410942237, 0.0198130534651379, 0.192904154459635, 1.31964428405416,
0.39110403527359, 0.332270985081281, 0.478428447791025, 0.0124228421288232,
0.108003384278466, 0.44015597643793, 0.0678367833438553, 0.340908625770017,
0.100580311380327, 0.241738101143046, 0.250514230575028, 0.0783424836236563,
0.359960379064981, 0.342748967916672, 0.430753882558053, 0.417701784817779,
0.18488046588997, 0.100427665458216, 0.431041552050316, 0.331518596024219,
0.17139143201833, 0.669277467426467, 0.140747481646637, 0.72625752274651,
1.07259633924935, 0.185943118296564, 0.19820588392516, 0.325798601941002,
0.0699555268511176, 0.103149285384764, 0.368645422767504, 0.258062588036958,
0.0298913594374683, 0.36939221892146, 0.0824214177800864, 0.523995615513929,
1.61093758162402, 0.143710710729162, 0.910129771215431, 0.378943805090106,
0.271122748654481, 0.279002163557281, 0.594921801513742, 0.770823875400138,
0.969962422990418, 0.0951969947976371, 0.129595590685436, 0.0173518159426749,
0.117290165961365, 0.521747115103197, 0.271511935442642, 0.919747929357203,
0.128731187824502, 0.336085485590899, 0.272838062945121, 0.0197537352020542,
0.761284491005254, 0.268056970386112, 0.527898693352688, 0.411263835931219,
0.448548006190264, 0.700125743282549, 0.345032382432028, 0.149150631080071,
0.347869504850014, 0.0869427483218412, 0.227076369337738, 0.0879127546757239,
0.148868551484028, 0.0702022976862887, 0.0441904686409075, 0.116296115132973
]
The likelihood of the data given the rate parameter $\lambda$ is as follows.
\begin{align} Pr(X|\lambda) & = \prod_{i=1}^{n=100} \lambda \rm{exp}(-\lambda x_i) \cr & = \lambda^n \rm{exp}(-\lambda n \bar{x}) \end{align}
where $\bar{x}$ is the sample mean of $X$.
As described in Wikipedia (which is remarkably handy for distributions), if we assume a Gamma($\alpha$,$\beta$) prior on the rate coefficient, the posterior distribution of $\lambda$ is Gamma($\alpha+n$,$\beta+n \bar{x}$), the conjugate prior for an exponential distribution.
\begin{align} Pr(\lambda|X,\alpha, \beta) & \propto Pr(X|\lambda) \times Pr(\lambda| \alpha, \beta) \cr & = \lambda^n \rm{exp}(-\lambda n \bar{x}) \times \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} \rm{exp}(-\lambda \beta) \cr & = \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha+n-1} \rm{exp}(-\lambda (\beta + n \bar{x})) \end{align}
The marginal likelihood of this model can be calculated by integrating $Pr(X|\lambda) \times Pr(\lambda| \alpha, \beta)$ over $\lambda$.
\begin{align} \int_{\lambda=0}^{\infty}\frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha+n-1} \rm{exp}(-\lambda (\beta + n \bar{x})) \; d\lambda & = \frac{\beta^\alpha}{\Gamma(\alpha)} \int_{0}^{\infty}\lambda^{\alpha+n-1} exp(-\lambda (\beta + n \bar{x})) \; d\lambda \cr & = \frac{\beta^\alpha}{\Gamma(\alpha)} \frac{\Gamma(\alpha+n)}{(\beta+ n \bar{x})^{a+n}} \end{align}
The log marginal likelihood in Julia can be calculated as follows.
function lml(x::Vector{Float64},alpha::Float64,beta::Float64)
mux=mean(x)
n=convert(Float64,length(x))
alpha*log(beta)-(alpha+n)*log(beta+n*mux)+lgamma(alpha+n)-lgamma(alpha)
end
For $\alpha=1$ and $\beta=1$, the log marginal likelihood for these data is around 3.6.
alpha=1.0
beta=1.0
lml(x,alpha,beta)
In many cases, however, we don't have an analytical solution to the posterior distribution or the marginal likelihood.
Tempering
Several approaches to calculating marginal likelhoods are based on the idea of tempering, in which we consider running MCMC at a range of different (inverse) 'temperatures', obtaining by raising likelihood to a power between 0 and 1; when the power is 0, we sample from the prior, while when the power is 1, we sample from the posterior. While we use the tempered likelihood to determine acceptance probabilities in the MCMC, we will use samples of the untempered likelihood to compute the marginal likelihood.
function met_temper(x::Vector{Float64},lambda0::Float64,alpha::Float64,beta::Float64,sigma::Float64,temp::Float64,niters::Int64)
lambdavec=vec(Array(Float64,niters))
llikvec=vec(Array(Float64,niters))
lliktvec=vec(Array(Float64,niters))
lpriorvec=vec(Array(Float64,niters))
lpostvec=vec(Array(Float64,niters))
lambda=lambda0
llik=sum(logpdf(Exponential(1.0/lambda),x))
llikt=temp*llik
lprior=logpdf(Gamma(alpha,beta),lambda)
lpost=llik+lprior
for i in 1:niters
lambdas = lambda + rand(Normal(0,sigma),1)[1]
if lambdas>0
lliks=sum(logpdf(Exponential(1/lambdas),x))
llikst=temp*lliks
lpriors=logpdf(Gamma(alpha,beta),lambdas)
lposts=llikst+lpriors
A=exp(lposts-lpost)
if rand()<A
lambda,llik,llikt,lprior,lpost=lambdas,lliks,llikst,lpriors,lposts
end
end
lambdavec[i],llikvec[i],lliktvec[i],lpriorvec[i],lpostvec[i]=lambda,llik,llikt,lprior,lpost
end
[lambdavec llikvec lliktvec lpriorvec lpostvec]
end
I first run the chain at a range of temperatures.
tempvec=linspace(0,1.0,101).^5.0 # Range of temperatures
numtemp=length(tempvec)
pplist=vec(Array(Matrix{Float64},numtemp))
burnin=1000
niters=10000
lambda0=1.0 # Initial condition for lambda
sigma=1.0 # Dispersion for the random walk sampler
for i in 1:numtemp
out=met_temper(x,lambda0,alpha,beta,sigma,tempvec[i],niters+burnin)
pplist[i]=out
end
Power posteriors
The power posterior approach is based on integrating the expectation of the log likelihood (for a given chain) across the inverse temperatures (see Friel and Pettit (2008) and Lartillot and Philippe (2006). Friel and Pettitt used a trapezoidal scheme to integrate the likelihood, while a revised scheme used an improved trapezoidal scheme, which also requires the variance of the log likelihood.
I first extract the mean and the variance of the log likelihoods.
ell=zeros(numtemp)
vll=zeros(numtemp)
for i in 1:numtemp
out=pplist[i]
ell[i]=mean(out[burnin:niters,2])
vll[i]=var(out[burnin:niters,2])
end
The following uses this information to compute the log marginal likelihood, using the modified trapezoidal scheme.
function ppml(ell::Vector{Float64},vll::Vector{Float64},tempvec::Vector{Float64})
N=length(ell)
res=0.0
for i in 1:(N-1)
wts = tempvec[i+1]-tempvec[i]
res = res+wts*((ell[i+1]+ell[i])/2.0)-((wts^2)/12.0)*(vll[i+1]-vll[i])
end
res
end
This function calculates the bounds on the marginal likelihood, based on the error from integration.
function boundml(ell::Vector{Float64},tempvec::Vector{Float64})
tempdiff = tempvec[2:length(tempvec)] .- tempvec[1:(length(tempvec)-1)]
ub = sum(tempdiff .* ell[2:length(ell)])
lb = sum(tempdiff .* ell[1:(length(ell)-1)])
[lb ub]
end
Now I can compute the log marginal likelihood from the series of tempered chains, and compare with the analytical result.
ppml(ell,vll,tempvec)
boundml(ell,tempvec)
Stepping stone
The stepping stone approach also employs tempered distributions, but rather than integrating the expectation of the log likelihood, it employs importance sampling between adjacent tempered chains to calculate a series of normalising constants. The product of these normalising constants gives an estimate of the marginal likelihood. Here, I estimate the log ratio of marginal likelihoods.
lrss=0.0
for i in 2:numtemp
tempdiff = tempvec[i]-tempvec[i-1]
logmaxll = maximum(pplist[i][(burnin+1):(burnin+niters),2])
oldll = pplist[i-1][(burnin+1):(burnin+niters),2]
lrss = lrss+tempdiff*logmaxll
lrss = lrss+log((1/length(oldll))*sum(exp(tempdiff*(oldll-logmaxll))))
end
lrss
As you can see, the code isn't that different between R and Julia (although the ability to do multiple assignment in a single line in Julia makes for more compact code). However, where they differ is in speed. The bottleneck is in running the MCMC across multiple temperatures. Firstly, let's calculate the time elapsed for a single run in Julia.
@elapsed out=met_temper(x,lambda0,alpha,beta,sigma,1.0,100000)
In R, the equivalent code runs in about 2.7 seconds, or about ten times slower.
No comments:
Post a Comment