# example 2 - a "non-parametric" regression
# with a locally quadratic prior
# source("example2.s")

 # make a fake dataset

n_50
sites_1:n
y_sin(sites/25*pi)+rnorm(n,sd=.2)
plot(sites,y)

 # now set some MCMC parameters

iter_600
burn_200
ay_1; by_.005; ax_1; bx_.005

 # make the matrix W

W_diag(rep(1,n)) 
for(i in 1:(n-1)){
   W[i,i+1]_-2/3
   W[i+1,i]_-2/3
}
for(i in 1:(n-2)){
   W[i,i+2]_1/6
   W[i+2,i]_1/6
}
for(i in 1:2) W[i,i]_-sum(W[i,-i])
for(i in (n-1):n) W[i,i]_-sum(W[i,-i])

 # arrays to hold the MCMC realizations

xout_matrix(NA,ncol=n,nrow=iter)
lamy_rep(NA,n)
lamx_rep(NA,n)

 # give initial values to parameters

xout[1,]_y
xout[1,]_lowess(sites,y)$y  # or smooth version of y
lamy[1]_10
lamx[1]_2

for(k in 2:iter){
  rss_sum((y-xout[k-1,])^2)
  lamy[k]_rgamma(1,ay+n/2)/(by+.5*rss)

  xss_t(xout[k-1,])%*%W%*%xout[k-1,]
  lamx[k]_rgamma(1,ax+n/2)/(bx+.5*xss)

  prec_diag(rep(lamy[k],n))+lamx[k]*W
  vmat_solve(prec)
  vmat_.5*(vmat + t(vmat))
  mean_as.vector(lamy[k]*vmat%*%y)
  xout[k,]_rmultnorm(1,mean,vmat)
}

 # now get posterior means...

xhat_apply(xout[-(1:burn),],2,mean)
lamyhat_mean(lamy[-(1:burn)])
lamxhat_mean(lamx[-(1:burn)])

 # plot the solution

par(mfrow=c(2,2),oma=c(0,0,0,0), mar=c(5,5,2,2))
plot(sites,y)
lines(sites,xhat)
xbounds_apply(xout[-(1:burn),],2,quantile,c(.1,.9))
matlines(sites,t(xbounds),lty=2)
mtext("posterior mean and 80% CI for x",side=3,line=1,cex=.8)

plot(lamy,type="l")
mtext("lamy",side=3,line=1,cex=.8)
plot(lamx,type="l")
mtext("lamx",side=3,line=1,cex=.8)
plot(xout[,40],type="l")
mtext("x[40]",side=3,line=1,cex=.8)