1 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)\).

2 Non-informative prior

The prior for parameters \((\alpha,\beta,\theta,\sigma^2)\) are as follows: \[ p(\alpha,\beta,\theta,\sigma^2) \propto \frac{1}{\sigma^2}, \] which represents a non-informative prior. Therefore, the joint posterior is \[ p(\alpha,\beta,\theta,\sigma^2|\mbox{data}) \propto \frac{1}{\sigma^2} \prod_{i=1}^n \left(\frac{1}{\sigma^2}\right)^{1/2}\exp\left\{-\frac{1}{2\sigma^2} \left(y_i-\alpha-\beta z_i(\theta)\right)^2 \right\}, \] where \(z_i(\theta) = x_i/(\theta+x_i)\) is explictly a function of \(\theta\). Posterior inference is approximated by a Markov chain Monte Carlo scheme that cycles through the full conditions of the four parameters, \(\alpha\), \(\beta\), \(\theta\) and \(\sigma^2\).

3 Full conditionals of \(\alpha\)

  1. 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\).

4 Full conditional of \(\beta\)

  1. 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] \]

5 Full conditional of \(\sigma^2\)

  1. 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
  2. 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.

6 Posterior predictive analysis

  1. 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=