Differences between revisions 2 and 3
Revision 2 as of 2009-10-26 17:21:22
Size: 265
Editor: mna75-13-88-175-67-73
Comment:
Revision 3 as of 2009-10-28 13:32:06
Size: 6151
Editor: mna75-13-88-175-67-73
Comment:
Deletions are marked like this. Additions are marked like this.
Line 10: Line 10:


== Introduction ==

Bayesian model selection is performed on the basis of the posterior model probabilities <<latex($p(k|y)$)>> which can be evaluate in this
case by considering the models parameters as unknown and updating at each iteration of the MCMC knowing the model index k.
To be more precise, the dimensionality of the parameter space represents an unknown parameter.


Let <<latex($\theta^k$)>> be the parameters vector associated with the model of order k.
There's three kind of possible move which can be described.

In steps 1 and 2, we go from a model k to an other one of different length k' with the probability <<latex($J(k'/k)$)>>.
The dimension matching requirement (Green, 1995) set dim(k)+dim(q(u))=dim(k')+dim(q(u')) where the function q(.) is the proposal
and u an additional random vector.
In the literature, we call a "birth move" the fact that dim(k')>dim(k), otherwise this is a "death move".

In the last step, we stay within a same model. In that case, a standard MH can be run.



== Algorithm ==

I. choose initial values

II. for i=1,N

{{{#!latex
$J(k'|k)\propto exp(\lambda)$
}}}

i. Birth move: if <<latex($dim(\theta^{(k')})>dim(\theta^{(k)})$)>>

go from k to k'=k+n by drawing n values from a random vector u of length n

{{{#!latex
\begin{eqnarray*}
\theta^{(k')}_u \sim p(\theta^{(k')}_u|k') \\
\theta^{(k')}=
\left (
   \begin{array}{c}
      \theta^{(k)} \\
      \theta^{(k')}_u \\
   \end{array}
   \right )
\end{eqnarray*}
}}}

if <<latex($x\sim U(0,1)<=\alpha^{birth}$)>>
       k=k'
end if

ii. Death move: else if <<latex($dim(\theta^{(k')})<dim(\theta^{(k)})$)>>

deterministic step to go from k to k'=k-n

<<latex($\theta^{k'}=\theta^{k}_{(1...k-n)}$)>>

if <<latex($x\sim U(0,1)<=\alpha^{death}$)>>
       k=k'
end if

iii. else <<latex($\theta^{k'}=\theta^k+N(0,\sigma)$)>>

                 end if

end


footnote:
The proposal distribution for changes in model order has the following expression:

{{{#!latex
$J(k'|k)\propto exp(\lambda)\frac{\lambda}{k!}$
}}}

or
{{{#!latex
$$J(k'|k)\propto \begin{cases}
       exp(-\lambda)|k'-k| & \text{if } k'\not=k \\
       0 & \text{if } k'=k
       \end{cases}
$$
}}}
The symmetry of the distribution leads the ratio of model prior densities to be 1

== Acceptation Probabilities ==
we assume here that the proposal distribution is the parameters prior density.

Consider

{{{#!latex
\begin{eqnarray*}
\theta^{'(i)}\sim(\theta_i|\theta_{-i},k,y)\propto
\begin{cases} p(\theta^{'(k')}|y,k',\theta^{(k)}) & \text{if } i=k'\\
p(\theta^{'(i)}) & \text{if } i\not= k'
\end{cases}
\end{eqnarray*}
}}}


where the proposal distribution is <<latex($ p(\theta^{'(k')}|y,k',\theta^{(k)})$)>> and <<latex($p(\theta^{'(i)})$)>>
is called "pseudo-prior".


Under the assumption of a symmetric random-walk proposal distribution, we can derive the acceptance ratio.


{{{#!latex
\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,\underbrace{\frac{p(k',\theta^{(k')}|y,\sigma_e)}{p(k,\theta^{(k)}|y,\sigma_e)}}_{target proba}
\underbrace{\frac{J(k|k')}{J(k'|k)}}_{model prior}
\underbrace{\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}}_{proposal ratio}
\right)
\end{eqnarray*}
}}}


The Jacobian cancel since the transformation let unchanged the scale of parameters vectors.


{{{#!latex
\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
\frac{p(k')}{p(k)}
\frac{p(\theta^{(k')}|k')}{p(\theta^{(k)}|k)}
\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}
\right)
\end{eqnarray*}
}}}

using the definition of joint posterior distribution:


{{{#!latex
\begin{eqnarray*}
p(k,\theta^{(k)},\sigma_e|y)\propto
\underbrace{p(y|k,\theta^{(k)},\sigma_e)}_{Likelihood}
\underbrace{p(k)p(\theta^{(k)}|k)p(\sigma_e)}_{priors}
\end{eqnarray*}
}}}

Precisely the full posterior distribution for the composite model space can be expressed as


{{{#!latex
\begin{eqnarray*}
p(k,\theta|y)\propto p(y|k,\theta)p(\theta_k|k)p(\theta_{-k}|\theta_k,k)p(k)
\end{eqnarray*}
}}}

It depends only upon the corresponding parameters <<latex($\theta_k$)>>, for a particular k.
<<latex($\theta_{-k}$)>> represents all parameters not used by model k. Hence, <<latex($p(\theta_{-k}|\theta_k,k)$)>> is this "pseudo-prior"
which is not used by k.


Assuming p(k) follows a uniform distribution, the ratio <<latex($\frac{p(k')}{p(k)}=1$)>> and
the common parameters prior in the ratio prior densities cancel. If we consider the case where <<latex($p(\theta^{(k)}|k)$)>> is iid gaussian,
the priors for the common model terms cancel leaving for a birth move:


{{{#!latex
\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
p_{\theta}(\theta_u^{(k')}|k')
%\frac{p(\theta^{(k')}|\sigma_e)_\rho}{p(\theta^{(k)}|\sigma_e)\rho}
\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}
\right)
\end{eqnarray*}
}}}

 Knowing that <<latex($q(\theta_u^{(k')}|\sigma_e)=p(\theta^{(k')}_u|\sigma_e)$)>> and the fact that
the reverse transformation is deterministic (ie <<latex($\theta^{(k)}_u$)>> is empty and <<latex($q(\theta^{(k)}_u)=1$))>>.
 We have finally:


{{{#!latex
\begin{eqnarray*}
\alpha ^{birth}\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
\right)
\end{eqnarray*}
}}}

and


{{{#!latex
\begin{eqnarray*}
\alpha ^{death}\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k,\theta^{(k)},\sigma_e)}{p(y|k',\theta^{(k')},\sigma_e)}
\right)
\end{eqnarray*}
}}}

Reversible Jump MCMC routines in Dynare

  • Functions

  • Structures

  • Notes on the algorithmrjmcmc

Introduction

Bayesian model selection is performed on the basis of the posterior model probabilities $p(k|y)$ which can be evaluate in this case by considering the models parameters as unknown and updating at each iteration of the MCMC knowing the model index k. To be more precise, the dimensionality of the parameter space represents an unknown parameter.

Let $\theta^k$ be the parameters vector associated with the model of order k. There's three kind of possible move which can be described.

In steps 1 and 2, we go from a model k to an other one of different length k' with the probability $J(k'/k)$. The dimension matching requirement (Green, 1995) set dim(k)+dim(q(u))=dim(k')+dim(q(u')) where the function q(.) is the proposal and u an additional random vector. In the literature, we call a "birth move" the fact that dim(k')>dim(k), otherwise this is a "death move".

In the last step, we stay within a same model. In that case, a standard MH can be run.

Algorithm

I. choose initial values

II. for i=1,N

$J(k'|k)\propto exp(\lambda)$

i. Birth move: if $dim(\theta^{(k')})>dim(\theta^{(k)})$

go from k to k'=k+n by drawing n values from a random vector u of length n

\begin{eqnarray*}
\theta^{(k')}_u \sim p(\theta^{(k')}_u|k') \\
\theta^{(k')}=
\left (
   \begin{array}{c}
      \theta^{(k)}  \\
      \theta^{(k')}_u   \\
   \end{array}
   \right )
\end{eqnarray*}

if $x\sim U(0,1)<=\alpha^{birth}$

  • k=k'

end if

ii. Death move: else if $dim(\theta^{(k')})<dim(\theta^{(k)})$

deterministic step to go from k to k'=k-n

$\theta^{k'}=\theta^{k}_{(1...k-n)}$

if $x\sim U(0,1)<=\alpha^{death}$

  • k=k'

end if

iii. else $\theta^{k'}=\theta^k+N(0,\sigma)$

  • end if

end

footnote: The proposal distribution for changes in model order has the following expression:

$J(k'|k)\propto exp(\lambda)\frac{\lambda}{k!}$

or

latex error! exitcode was 1 (signal 0), transscript follows:

This is pdfTeX, Version 3.14159265-2.6-1.40.17 (TeX Live 2016/Debian) (preloaded format=latex)
entering extended mode
(./latex_1d4470da2edd2c07fe2b2d127bf8617649383edb_p.tex
LaTeX2e <2017/01/01> patch level 3
Babel <3.9r> and hyphenation patterns for 3 language(s) loaded.
(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
(/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo))
(/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty
(/usr/share/texlive/texmf-dist/tex/latex/base/utf8.def
(/usr/share/texlive/texmf-dist/tex/latex/base/t1enc.dfu)
(/usr/share/texlive/texmf-dist/tex/latex/base/ot1enc.dfu)
(/usr/share/texlive/texmf-dist/tex/latex/base/omsenc.dfu)))
No file latex_1d4470da2edd2c07fe2b2d127bf8617649383edb_p.aux.
! Misplaced alignment tab character &.
l.8        exp(-\lambda)|k'-k| &
                                 \text{if } k'\not=k \\
! Undefined control sequence.
l.8        exp(-\lambda)|k'-k| & \text
                                      {if } k'\not=k \\
! Misplaced alignment tab character &.
l.9        0 &
               \text{if } k'=k
! Undefined control sequence.
l.9        0 & \text
                    {if } k'=k
[1] (./latex_1d4470da2edd2c07fe2b2d127bf8617649383edb_p.aux) )
(see the transcript file for additional information)
Output written on latex_1d4470da2edd2c07fe2b2d127bf8617649383edb_p.dvi (1 page,
 500 bytes).
Transcript written on latex_1d4470da2edd2c07fe2b2d127bf8617649383edb_p.log.

The symmetry of the distribution leads the ratio of model prior densities to be 1

Acceptation Probabilities

we assume here that the proposal distribution is the parameters prior density.

Consider

latex error! exitcode was 1 (signal 0), transscript follows:

This is pdfTeX, Version 3.14159265-2.6-1.40.17 (TeX Live 2016/Debian) (preloaded format=latex)
entering extended mode
(./latex_65910c034e45627534ec6937ef8e587fe68d960e_p.tex
LaTeX2e <2017/01/01> patch level 3
Babel <3.9r> and hyphenation patterns for 3 language(s) loaded.
(/usr/share/texlive/texmf-dist/tex/latex/base/article.cls
Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
(/usr/share/texlive/texmf-dist/tex/latex/base/size12.clo))
(/usr/share/texlive/texmf-dist/tex/latex/base/inputenc.sty
(/usr/share/texlive/texmf-dist/tex/latex/base/utf8.def
(/usr/share/texlive/texmf-dist/tex/latex/base/t1enc.dfu)
(/usr/share/texlive/texmf-dist/tex/latex/base/ot1enc.dfu)
(/usr/share/texlive/texmf-dist/tex/latex/base/omsenc.dfu)))
No file latex_65910c034e45627534ec6937ef8e587fe68d960e_p.aux.
! Extra }, or forgotten \endgroup.
<template> }
            $\@eqnsel \endtemplate 
l.9 ...ases} p(\theta^{'(k')}|y,k',\theta^{(k)}) &
                                                   \text{if } i=k'\\
! Missing \endgroup inserted.
<inserted text> 
                \endgroup 
l.9 ...ases} p(\theta^{'(k')}|y,k',\theta^{(k)}) &
                                                   \text{if } i=k'\\
! Missing } inserted.
<inserted text> 
                }
l.9 ...ases} p(\theta^{'(k')}|y,k',\theta^{(k)}) &
                                                   \text{if } i=k'\\
! Undefined control sequence.
l.9 ...p(\theta^{'(k')}|y,k',\theta^{(k)}) & \text
                                                  {if } i=k'\\
! Undefined control sequence.
l.10 p(\theta^{'(i)}) & \text
                             {if } i\not= k'

! LaTeX Error: \begin{eqnarray*} on input line 7 ended by \end{cases}.

See the LaTeX manual or LaTeX Companion for explanation.
Type  H <return>  for immediate help.
 ...                                              
                                                  
l.11 \end{cases}
                
! Missing } inserted.
<inserted text> 
                }
l.11 \end{cases}
                
! Missing $ inserted.
<inserted text> 
                $
l.11 \end{cases}
                
! Missing } inserted.
<inserted text> 
                }
l.11 \end{cases}
                
! Missing \cr inserted.
<inserted text> 
                \cr 
l.11 \end{cases}
                
! Missing { inserted.
<inserted text> 
                {
l.11 \end{cases}
                
! Missing $ inserted.
<inserted text> 
                $
l.11 \end{cases}
                
! Missing $$ inserted.
<to be read again> 
                   \endgroup 
l.11 \end{cases}
                
! Misplaced alignment tab character &.
\reserved@a ->&
                &
l.12 \end{eqnarray*}
                    
! Misplaced alignment tab character &.
\reserved@a ->& &
                 
l.12 \end{eqnarray*}
                    
! Misplaced \cr.
\@@eqncr ...l \@eqnswtrue \global \@eqcnt \z@ \cr 
                                                  
l.12 \end{eqnarray*}
                    
! Too many }'s.
\endeqnarray ->\@@eqncr \egroup 
                                \global \advance \c@equation \m@ne $$\@ignor...
l.12 \end{eqnarray*}
                    

! LaTeX Error: \begin{document} ended by \end{eqnarray*}.

See the LaTeX manual or LaTeX Companion for explanation.
Type  H <return>  for immediate help.
 ...                                              
                                                  
l.12 \end{eqnarray*}
                    
! Missing $ inserted.
<inserted text> 
                $
l.12 \end{eqnarray*}
                    
! Display math should end with $$.
<to be read again> 
                   \endgroup 
l.12 \end{eqnarray*}
                    
! Extra \endgroup.
<recently read> \endgroup 
                          
l.12 \end{eqnarray*}
                    
[1] (./latex_65910c034e45627534ec6937ef8e587fe68d960e_p.aux) )
(see the transcript file for additional information)
Output written on latex_65910c034e45627534ec6937ef8e587fe68d960e_p.dvi (1 page,
 824 bytes).
Transcript written on latex_65910c034e45627534ec6937ef8e587fe68d960e_p.log.

where the proposal distribution is $ p(\theta^{'(k')}|y,k',\theta^{(k)})$ and $p(\theta^{'(i)})$ is called "pseudo-prior".

Under the assumption of a symmetric random-walk proposal distribution, we can derive the acceptance ratio.

\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,\underbrace{\frac{p(k',\theta^{(k')}|y,\sigma_e)}{p(k,\theta^{(k)}|y,\sigma_e)}}_{target  proba}
\underbrace{\frac{J(k|k')}{J(k'|k)}}_{model prior}
\underbrace{\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}}_{proposal ratio}
\right)
\end{eqnarray*}

The Jacobian cancel since the transformation let unchanged the scale of parameters vectors.

\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
\frac{p(k')}{p(k)}
\frac{p(\theta^{(k')}|k')}{p(\theta^{(k)}|k)}
\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}
\right)
\end{eqnarray*}

using the definition of joint posterior distribution:

\begin{eqnarray*}
p(k,\theta^{(k)},\sigma_e|y)\propto 
\underbrace{p(y|k,\theta^{(k)},\sigma_e)}_{Likelihood}
\underbrace{p(k)p(\theta^{(k)}|k)p(\sigma_e)}_{priors}
\end{eqnarray*}

Precisely the full posterior distribution for the composite model space can be expressed as

\begin{eqnarray*}
p(k,\theta|y)\propto p(y|k,\theta)p(\theta_k|k)p(\theta_{-k}|\theta_k,k)p(k)
\end{eqnarray*}

It depends only upon the corresponding parameters $\theta_k$, for a particular k. $\theta_{-k}$ represents all parameters not used by model k. Hence, $p(\theta_{-k}|\theta_k,k)$ is this "pseudo-prior" which is not used by k.

Assuming p(k) follows a uniform distribution, the ratio $\frac{p(k')}{p(k)}=1$ and the common parameters prior in the ratio prior densities cancel. If we consider the case where $p(\theta^{(k)}|k)$ is iid gaussian, the priors for the common model terms cancel leaving for a birth move:

\begin{eqnarray*}
\alpha\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
p_{\theta}(\theta_u^{(k')}|k')
%\frac{p(\theta^{(k')}|\sigma_e)_\rho}{p(\theta^{(k)}|\sigma_e)\rho}
\frac{q(\theta^{(k')}|k')}{q(\theta^{(k)}|k)}
\right)
\end{eqnarray*}

  • Knowing that $q(\theta_u^{(k')}|\sigma_e)=p(\theta^{(k')}_u|\sigma_e)$ and the fact that

the reverse transformation is deterministic (ie $\theta^{(k)}_u$ is empty and $q(\theta^{(k)}_u)=1$).

  • We have finally:

\begin{eqnarray*}
\alpha ^{birth}\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k',\theta^{(k')},\sigma_e)}{p(y|k,\theta^{(k)},\sigma_e)}
\right)
\end{eqnarray*}

and

\begin{eqnarray*}
\alpha ^{death}\left((k,\theta^{(k)}\to (k',\theta^{(k')})\right)=
min\left(1,
\frac{p(y|k,\theta^{(k)},\sigma_e)}{p(y|k',\theta^{(k')},\sigma_e)}
\right)
\end{eqnarray*}

DynareWiki: ReversibleJumpMCMC (last edited 2018-10-17 15:56:59 by StéphaneAdjemian)