# example 1 - a "non-parametric" regression
# source("example1.s")

 # make a fake dataset

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

 # now set some MCMC parameters

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

 # make the matrix W

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

 # 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]_1
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_sum(diff(xout[k-1,])^2)
  # or, less efficiently: 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...

burn_200
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)