Model
The following example was taken from Chapter 6 (Metropolis-Hastings
algorithms) of our book Gamerman and Lopes (2006) Markov Chain
Monte Carlo: Stochastic Simulation for Bayesian Inference.
The variable \(x_i\) is used to
explain \(y_i\) on a nonlinear Gaussian
regression set up: \[
y_i|x_i,\alpha,\beta,\theta,\sigma^2 \sim N(y_i,\alpha+\beta
x_i/(\theta+x_i),\sigma^2).
\] The observed data is data=\(\{x,y\}\), where \(x=(0.02,0.02,0.06,0.06,0.11,0.11,0.22,0.22,0.56,0.56,1.10,1.10)\)
and \(y=(76,47,97,107,123,139,159,152,191,201,207,200)\).
Full conditionals of
\(\alpha\)
- Show that the full conditional distribution of \(\alpha\) is as follows: \[
\alpha|\beta,\theta,\sigma^2,\mbox{data} \sim N\left[
{\bar y}-\beta {\bar z}(\theta),\sigma^2/n
\right],
\] where \({\bar y}=\sum_{i=1}^n
y_i/n\) and \({\bar
z}(\theta)=\sum_{i=1}^n z_i(\theta)/n\).
Full conditional of
\(\beta\)
- Show that the full conditional distribution of \(\beta\) is as follows: \[
\beta|\alpha,\theta,\sigma^2,\mbox{data} \sim
N\left[\frac{\sum_{i=1}^n (y_i-\alpha)z_i(\theta)}{\sum_{i=1}^n
z_i^2(\theta)},
\frac{\sigma^2}{\sum_{i=1}^n z_i^2(\theta)}\right]
\]
Full conditional of
\(\sigma^2\)
- Show that the full conditional distribution of \(\sigma^2\) is as follows: \[
\sigma^2|\alpha,\beta,\theta,\mbox{data} \sim
IG\left[\frac{n-3}{2},\frac{
\sum_{i=1}^n (y_i-\alpha-\beta z_i(\theta))^2.
}{2}\right]
\] # Implementing the MCMC algorithm
- Implement the MCMC scheme and verify the trace plots for mixing and
convergence of the chains. Also, plot the autocorrelation functions of
the sequences of simulated draws. Finally, plot the histograms of the
marginal posterior distributions for each one of the four parameters.
The parameter \(\theta\) is the only
one whose full conditional distribution does not belong to any known
family of distributions, so we will use a Metropolis-Hastings within our
MCMC scheme to sample from \(\theta\)
given \((\alpha,\beta,\sigma^2)\) and
the data. For the random-walk Metropolis proposal we will use \(\theta^* \sim N(\theta^{old},\tau^2)\),
where \(\tau=0.01\) is the tuning
variance of the proposal. Aditionally, we will start our MCMC scheme
with \(\alpha^{(0)}=50\), \(\beta^{(0)}=170\), \(\theta^{(0)}=0.13\) and \(\sigma^2=126\). Run your MCMC for a burn-in
period of \(M_0=10000\) draws and then
keep the next \(M=10000\) draws from
posterior summaries and inference.
Posterior predictive
analysis
- For \(x_{new} \in
\{0.02,0.03,\ldots,1.10\}\), obtain MCMC approximations to $
\[
E(y_{new}|x_{new},\mbox{data}) = \int
E(y_{new}|x_{new},\alpha,\beta,\theta)p(\alpha,\beta,\theta|\mbox{data})d\alpha
d\beta d\theta,
\] where \[
E_{MC}(y_{new}|x_{new},\mbox{data}) = \frac{1}{M} \sum_{i=1}^M
E(y_{new}|x_{new},\alpha^{(i)},\beta^{(i)},\theta^{(i)}) =
\frac{1}{M} \sum_{i=1}^M \left(\alpha^{(i)}+\beta^{(i)}
\frac{x_{new}}{\theta^{(i)}+x_{new}}\right),
\] where \((\alpha^{(i)},\beta^{(i)},\theta^{(i)})\)
are the MCMC-based posterior draws from \(p(\alpha,\beta,\theta|\mbox{data})\). Plot
the scatterplot of \(x_i\) vs \(y_i\) (observed data) and include in the
picture the pairs \((x_{new},E(y_{new}|x_{new},\mbox{data}))\).
Notice the nonlinear relationship.
LS0tCnRpdGxlOiAiQWR2YW5jZWQgQmF5ZXNpYW4gU3RhdGlzdGljYWwgTGVhcm5pbmciCnN1YnRpdGxlOiAiSFczOiBHaWJicyBhbmQgTWV0cm9wb2xpcyBzdGVwcyIKYXV0aG9yOiAiSGVkaWJlcnQgRnJlaXRhcyBMb3BlcyIKZGF0ZTogIjIvMjQvMjAyMyIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6CiAgICB0b2M6IHRydWUKICAgIHRvY19kZXB0aDogMgogICAgY29kZV9kb3dubG9hZDogeWVzCiAgICBudW1iZXJfc2VjdGlvbnM6IHRydWUKLS0tCiAgCiAgYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CmtuaXRyOjpvcHRzX2NodW5rJHNldChlY2hvID0gVFJVRSkKYGBgCgojIE1vZGVsClRoZSBmb2xsb3dpbmcgZXhhbXBsZSB3YXMgdGFrZW4gZnJvbSBDaGFwdGVyIDYgKE1ldHJvcG9saXMtSGFzdGluZ3MgYWxnb3JpdGhtcykgb2Ygb3VyIGJvb2sgR2FtZXJtYW4gYW5kIExvcGVzICgyMDA2KSAqKk1hcmtvdiBDaGFpbiBNb250ZSBDYXJsbzogU3RvY2hhc3RpYyBTaW11bGF0aW9uIGZvciBCYXllc2lhbiBJbmZlcmVuY2UqKi4gIAoKVGhlIHZhcmlhYmxlICR4X2kkIGlzIHVzZWQgdG8gZXhwbGFpbiAkeV9pJCBvbiBhIG5vbmxpbmVhciBHYXVzc2lhbiByZWdyZXNzaW9uIHNldCB1cDoKJCQKeV9pfHhfaSxcYWxwaGEsXGJldGEsXHRoZXRhLFxzaWdtYV4yIFxzaW0gTih5X2ksXGFscGhhK1xiZXRhIHhfaS8oXHRoZXRhK3hfaSksXHNpZ21hXjIpLgokJApUaGUgb2JzZXJ2ZWQgZGF0YSBpcyBkYXRhPSRce3gseVx9JCwgd2hlcmUgJHg9KDAuMDIsMC4wMiwwLjA2LDAuMDYsMC4xMSwwLjExLDAuMjIsMC4yMiwwLjU2LDAuNTYsMS4xMCwxLjEwKSQgYW5kICR5PSg3Niw0Nyw5NywxMDcsMTIzLDEzOSwxNTksMTUyLDE5MSwyMDEsMjA3LDIwMCkkLgoKIyBOb24taW5mb3JtYXRpdmUgcHJpb3IKVGhlIHByaW9yIGZvciBwYXJhbWV0ZXJzICQoXGFscGhhLFxiZXRhLFx0aGV0YSxcc2lnbWFeMikkIGFyZSBhcyBmb2xsb3dzOgokJApwKFxhbHBoYSxcYmV0YSxcdGhldGEsXHNpZ21hXjIpIFxwcm9wdG8gXGZyYWN7MX17XHNpZ21hXjJ9LAokJAp3aGljaCByZXByZXNlbnRzIGEgbm9uLWluZm9ybWF0aXZlIHByaW9yLiAgVGhlcmVmb3JlLCB0aGUgam9pbnQgcG9zdGVyaW9yIGlzCiQkCnAoXGFscGhhLFxiZXRhLFx0aGV0YSxcc2lnbWFeMnxcbWJveHtkYXRhfSkgXHByb3B0byBcZnJhY3sxfXtcc2lnbWFeMn0KXHByb2Rfe2k9MX1ebiBcbGVmdChcZnJhY3sxfXtcc2lnbWFeMn1ccmlnaHQpXnsxLzJ9XGV4cFxsZWZ0XHstXGZyYWN7MX17MlxzaWdtYV4yfSBcbGVmdCh5X2ktXGFscGhhLVxiZXRhIHpfaShcdGhldGEpXHJpZ2h0KV4yIFxyaWdodFx9LAokJAp3aGVyZSAkel9pKFx0aGV0YSkgPSB4X2kvKFx0aGV0YSt4X2kpJCBpcyBleHBsaWN0bHkgYSBmdW5jdGlvbiBvZiAkXHRoZXRhJC4gIFBvc3RlcmlvciBpbmZlcmVuY2UgaXMgYXBwcm94aW1hdGVkIGJ5IGEgTWFya292IGNoYWluIE1vbnRlIENhcmxvIHNjaGVtZSB0aGF0IGN5Y2xlcyB0aHJvdWdoIHRoZSBmdWxsIGNvbmRpdGlvbnMgb2YgdGhlIGZvdXIgcGFyYW1ldGVycywgJFxhbHBoYSQsICRcYmV0YSQsICRcdGhldGEkIGFuZCAkXHNpZ21hXjIkLgoKIyBGdWxsIGNvbmRpdGlvbmFscyBvZiAkXGFscGhhJAoKYSkgU2hvdyB0aGF0IHRoZSBmdWxsIGNvbmRpdGlvbmFsIGRpc3RyaWJ1dGlvbiBvZiAkXGFscGhhJCBpcyBhcyBmb2xsb3dzOgokJApcYWxwaGF8XGJldGEsXHRoZXRhLFxzaWdtYV4yLFxtYm94e2RhdGF9IFxzaW0gTlxsZWZ0Wwp7XGJhciB5fS1cYmV0YSB7XGJhciB6fShcdGhldGEpLFxzaWdtYV4yL24KXHJpZ2h0XSwKJCQKd2hlcmUgJHtcYmFyIHl9PVxzdW1fe2k9MX1ebiB5X2kvbiQgYW5kICR7XGJhciB6fShcdGhldGEpPVxzdW1fe2k9MX1ebiB6X2koXHRoZXRhKS9uJC4KCiMgRnVsbCBjb25kaXRpb25hbCBvZiAkXGJldGEkCmIpIFNob3cgdGhhdCB0aGUgZnVsbCBjb25kaXRpb25hbCBkaXN0cmlidXRpb24gb2YgJFxiZXRhJCBpcyBhcyBmb2xsb3dzOgokJApcYmV0YXxcYWxwaGEsXHRoZXRhLFxzaWdtYV4yLFxtYm94e2RhdGF9IFxzaW0gCk5cbGVmdFtcZnJhY3tcc3VtX3tpPTF9Xm4gKHlfaS1cYWxwaGEpel9pKFx0aGV0YSl9e1xzdW1fe2k9MX1ebiB6X2leMihcdGhldGEpfSwKXGZyYWN7XHNpZ21hXjJ9e1xzdW1fe2k9MX1ebiB6X2leMihcdGhldGEpfVxyaWdodF0KJCQKCiMgRnVsbCBjb25kaXRpb25hbCBvZiAkXHNpZ21hXjIkCmMpIFNob3cgdGhhdCB0aGUgZnVsbCBjb25kaXRpb25hbCBkaXN0cmlidXRpb24gb2YgJFxzaWdtYV4yJCBpcyBhcyBmb2xsb3dzOgokJApcc2lnbWFeMnxcYWxwaGEsXGJldGEsXHRoZXRhLFxtYm94e2RhdGF9IFxzaW0gSUdcbGVmdFtcZnJhY3tuLTN9ezJ9LFxmcmFjewpcc3VtX3tpPTF9Xm4gKHlfaS1cYWxwaGEtXGJldGEgel9pKFx0aGV0YSkpXjIuCn17Mn1ccmlnaHRdCiQkCiMgSW1wbGVtZW50aW5nIHRoZSBNQ01DIGFsZ29yaXRobQpkKSBJbXBsZW1lbnQgdGhlIE1DTUMgc2NoZW1lIGFuZCB2ZXJpZnkgdGhlIHRyYWNlIHBsb3RzIGZvciBtaXhpbmcgYW5kIGNvbnZlcmdlbmNlIG9mIHRoZSBjaGFpbnMuIEFsc28sIHBsb3QgdGhlIGF1dG9jb3JyZWxhdGlvbiBmdW5jdGlvbnMgb2YgdGhlIHNlcXVlbmNlcyBvZiBzaW11bGF0ZWQgZHJhd3MuICBGaW5hbGx5LCBwbG90IHRoZSBoaXN0b2dyYW1zIG9mIHRoZSBtYXJnaW5hbCBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9ucyBmb3IgZWFjaCBvbmUgb2YgdGhlIGZvdXIgcGFyYW1ldGVycy4gIFRoZSBwYXJhbWV0ZXIgJFx0aGV0YSQgaXMgdGhlIG9ubHkgb25lIHdob3NlIGZ1bGwgY29uZGl0aW9uYWwgZGlzdHJpYnV0aW9uIGRvZXMgbm90IGJlbG9uZyB0byBhbnkga25vd24gZmFtaWx5IG9mIGRpc3RyaWJ1dGlvbnMsIHNvIHdlIHdpbGwgdXNlIGEgTWV0cm9wb2xpcy1IYXN0aW5ncyB3aXRoaW4gb3VyIE1DTUMgc2NoZW1lIHRvIHNhbXBsZSBmcm9tICRcdGhldGEkIGdpdmVuICQoXGFscGhhLFxiZXRhLFxzaWdtYV4yKSQgYW5kIHRoZSBkYXRhLiAgRm9yIHRoZSByYW5kb20td2FsayBNZXRyb3BvbGlzIHByb3Bvc2FsIHdlIHdpbGwgdXNlICRcdGhldGFeKiBcc2ltIE4oXHRoZXRhXntvbGR9LFx0YXVeMikkLCB3aGVyZSAkXHRhdT0wLjAxJCBpcyB0aGUgdHVuaW5nIHZhcmlhbmNlIG9mIHRoZSBwcm9wb3NhbC4gIEFkaXRpb25hbGx5LCB3ZSB3aWxsIHN0YXJ0IG91ciBNQ01DIHNjaGVtZSB3aXRoICRcYWxwaGFeeygwKX09NTAkLCAkXGJldGFeeygwKX09MTcwJCwgCiRcdGhldGFeeygwKX09MC4xMyQgYW5kICRcc2lnbWFeMj0xMjYkLiAgUnVuIHlvdXIgTUNNQyBmb3IgYSBidXJuLWluIHBlcmlvZCBvZiAkTV8wPTEwMDAwJCBkcmF3cyBhbmQgdGhlbiBrZWVwIHRoZSBuZXh0ICRNPTEwMDAwJCBkcmF3cyBmcm9tIHBvc3RlcmlvciBzdW1tYXJpZXMgYW5kIGluZmVyZW5jZS4KCiMgUG9zdGVyaW9yIHByZWRpY3RpdmUgYW5hbHlzaXMgCmUpIEZvciAkeF97bmV3fSBcaW4gXHswLjAyLDAuMDMsXGxkb3RzLDEuMTBcfSQsIG9idGFpbiBNQ01DIGFwcHJveGltYXRpb25zIHRvICQKJCQKRSh5X3tuZXd9fHhfe25ld30sXG1ib3h7ZGF0YX0pID0gXGludCBFKHlfe25ld318eF97bmV3fSxcYWxwaGEsXGJldGEsXHRoZXRhKXAoXGFscGhhLFxiZXRhLFx0aGV0YXxcbWJveHtkYXRhfSlkXGFscGhhIGRcYmV0YSBkXHRoZXRhLAokJAp3aGVyZSAKJCQKRV97TUN9KHlfe25ld318eF97bmV3fSxcbWJveHtkYXRhfSkgPSBcZnJhY3sxfXtNfSBcc3VtX3tpPTF9Xk0gRSh5X3tuZXd9fHhfe25ld30sXGFscGhhXnsoaSl9LFxiZXRhXnsoaSl9LFx0aGV0YV57KGkpfSkgPSAKXGZyYWN7MX17TX0gXHN1bV97aT0xfV5NIFxsZWZ0KFxhbHBoYV57KGkpfStcYmV0YV57KGkpfSBcZnJhY3t4X3tuZXd9fXtcdGhldGFeeyhpKX0reF97bmV3fX1ccmlnaHQpLAokJAp3aGVyZSAkKFxhbHBoYV57KGkpfSxcYmV0YV57KGkpfSxcdGhldGFeeyhpKX0pJCBhcmUgdGhlIE1DTUMtYmFzZWQgcG9zdGVyaW9yIGRyYXdzIGZyb20gJHAoXGFscGhhLFxiZXRhLFx0aGV0YXxcbWJveHtkYXRhfSkkLiAgUGxvdCB0aGUgc2NhdHRlcnBsb3Qgb2YgJHhfaSQgdnMgJHlfaSQgKG9ic2VydmVkIGRhdGEpIGFuZCBpbmNsdWRlIGluIHRoZSBwaWN0dXJlIHRoZSBwYWlycyAkKHhfe25ld30sRSh5X3tuZXd9fHhfe25ld30sXG1ib3h7ZGF0YX0pKSQuICBOb3RpY2UgdGhlIG5vbmxpbmVhciByZWxhdGlvbnNoaXAuCgo=