# example 2 - a "non-parametric" regression
# with a Gaussian field prior
# source("example3.s")

 # make a fake dataset

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

 # make a function that gives covariance as a function
 # of distance

covpow_function(locs,pow=2,scale=5){
 # makes a covariance matrix using the power
 # exponential family C(d) = exp(-|d/scale|^pow)
 # locs is assumed to be a n x 2 matrix with each
 # row corresponding to locations
  # browser()
  d1_dist(locs)
  n_dim(locs)[1]
  mat_matrix(0,n,n)
  mat[lower.tri(mat)]_d1
  mat_mat+t(mat)
  cc_exp(-(mat/scale)^pow)
  return(cc)
}
svdinv_function(a,tol=1e-8){
 sva_svd(a)
 maxsv_max(sva$d)
 dd_ifelse(sva$d<tol*maxsv,0,1/sva$d)
 t(sva$u %*% (t(sva$v)*dd))
}


 # now set some MCMC parameters

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

 # make the matrix W

locs_cbind(sites,0)
Winv_covpow(locs,pow=2,scale=10)+diag(rep(.001,length(sites)))
W_svdinv(Winv)
# look at the full conditional of the 25th element:
# W[25,]
# [1]    0.1504459    0.5204976   -1.0747180   -0.9991295    1.2306029
# [6]    2.3462070   -0.2033997   -3.7420963   -2.6579320    3.5268332
#[11]    6.9253853    0.4001487  -10.1643600   -9.1589530    7.7229203
#[16]   20.6108293    5.8490640  -27.0699982  -33.3204781   14.0381100
#[21]   70.4320660   44.2085076 -104.2014056 -294.3764673  618.1652972
#[26] -294.3789562 -104.2049150   44.2076950   70.4367214   14.0443364
#[31]  -33.3217045  -27.0808292    5.8412274   20.6205487    7.7432864
#[36]   -9.1563116  -10.1931135    0.3717193    6.9452142    3.5872243
#[41]   -2.6374470   -3.8189291   -0.3005757    2.3888957    1.4280210
#[46]   -0.9153163   -1.3883893    0.1338830    1.0846310   -0.4629621
# sum(W[25,]) is .134.  This means there's shrinkage towards 0

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