### All required functions for using ESS ### Functions related to Wood and Chan algorithm of drawing samples ### Covariance matrix and design matrix (using basis function) are also defined! ### Required libraries: library(fields) library(FastGP) library(mcmcse) ### Define design matrix ### ### The basis functions ### # design matrix for monotone function estimation: psi_j=function(x,my_knot,delta_N) { N=length(my_knot) k=rep(0,N) i=max(which(my_knot<=x)) if(i==1) { k[1]=x-0.5*(x^2)/delta_N k[2]=x-my_knot[2]*x/delta_N+0.5*x^2/delta_N } if(i==2) { k[1]=delta_N/2 k[2]=delta_N/2+(x-my_knot[2])*(1+my_knot[2]/delta_N)-0.5*(x^2-my_knot[2]^2)/delta_N k[3]=(x-my_knot[2])*(1-my_knot[3]/delta_N)+0.5*(x^2-my_knot[2]^2)/delta_N } if(i==N) { k[1]=delta_N/2 k[2:(N-1)]=delta_N k[N]=delta_N/2 } if(i!=1 && i!=2 && i!=N) { k[1]=delta_N/2 k[2:(i-1)]=delta_N k[i]=delta_N/2+(x-my_knot[i])*(1+my_knot[i]/delta_N)-0.5*(x^2-my_knot[i]^2)/delta_N k[i+1]=(x-my_knot[i])*(1-my_knot[i+1]/delta_N)+0.5*(x^2-my_knot[i]^2)/delta_N } return(k) } ### Function to form design matrix: des.mat1=function(x,my_knot,delta_N){ n=length(x) N=length(my_knot)-1 # design matrix \Psi(n X N+1) X=matrix(0,n,N+1) for(l in 1:n){ X[l,1:(N+1)]=psi_j(x[l],my_knot,delta_N) } return(X) } #Given a \nu (smoothness parameter of Matern kernel) finding a value of # l (length-scale parameter) such that the correlation between the # maximum separation is some small value, say 0.05 #Matern kernel with smoothness nu and length-scale l: MK = function(x, y ,l, nu){ ifelse(abs(x-y)>0, (sqrt(2*nu)*abs(x-y)/l)^nu/(2^(nu-1)*gamma(nu))*besselK(x=abs(x-y)*sqrt(2*nu)/l, nu=nu), 1.0) } # function for uniroot using MK: fl=function(l,para){ #para[1]=x, para[2]=y and para[3]=nu of MK : Matern kernel function; #para[4]=pre-specified value of the correlation a=MK(para[1],para[2],l,para[3]) return(a-para[4]) } # function for estimating l for MK: l_est=function(nu,range,val){ # nu : smoothness; range : c(min, max) of the range of variable # val : pre-specified value of the correlation between the maximum separation para=c(range[1],range[2],nu,val) rl=uniroot(f=fl,interval = c(0.000001,100000),para) return(rl$root) } # Finding a value of l (length-scale parameter) such that the correlation between # the maximum separation is some small value, say 0.05 #Squared exponential kernel with length-scale l: SEK <- function(x,y,l){ ifelse(abs(x-y)>0, exp(- 0.5 * ((x-y)/l)^2), 1.0 + 10^(-8)) #exp(- 0.5 * ((x-y)/l)^2) } # function for uniroot using SEK: fl2=function(l,para){ #para[1]=x, para[2]=y, para[3]=pre-specified value of the correlation a=SEK(para[1],para[2],l) return(a-para[3]) } # function for estimating l for SEK: l_est2=function(range,val){ # range : c(min, max) of the range of variable # val : pre-specified value of the correlation between the maximum seperation para=c(range[1],range[2],val) rl=uniroot(f=fl2,interval = c(0.000001,100000),para) return(rl$root) } # Covariance matrix # MK covmat=function(knot,nu,l){ return(MK(rdist(knot),0,l,nu)) } # SEK covmat.sek=function(knot,l){ return(SEK(rdist(knot),0,l)) } # Order of the circulant matrix: # minimum value of g and m so that G can be embedded into C min_g=function(knot){ N=length(knot) g=ceiling(log(2*N,2)) #m=2^g and m>=2(n-1) : Wood & Chan notation; #since we are going upto n and not stopping at (n-1), the condition is modified! return("g" = g) } # forming the circulant matrix: circulant=function(x){ n = length(x) mat = matrix(0, n, n) for (j in 1:n) { mat[j, ] <- c(x[-(1:(n+1-j))], x[1:(n+1-j)]) } return(mat) } # Function for forming the vector of circulant matrix: circ_vec=function(knot,g,nu,l,tausq,mk=TRUE){ delta_N=1/(length(knot)-1) m=2**g cj=integer() for(j in 1:m){ if(j<=(m/2)) cj[j]=(j-1)*delta_N else cj[j]=(m-(j-1))*delta_N } if(mk == TRUE){ x=(tausq*MK(cj,0,l,nu)) }else{ x=(tausq*SEK(cj,0,l)) } return(x) } # Function for finding a g such that C is nnd: eig.eval=function(knot,g,nu,l,tausq,mk=TRUE){ vec=circ_vec(knot,g,nu,l,tausq,mk) C=circulant(vec) ev=min(eigen(C)$values) return(list("vec" = vec, "min.eig.val" = ev)) } # Function for finding a g such that C is nnd: # without forming the circulant matrix and without computing eigen values: C.eval=function(knot,g,nu,l,tausq,mk=TRUE){ vec=circ_vec(knot,g,nu,l,tausq,mk) val=fft(vec) # eigenvalues will be real as the circulant matrix formed by the # vector is by construction is symmetric! ev=min(Re(val)) return(list("vec" = vec, "min.eig.val" = ev)) } # Finding the smallest g such that C is nnd: nnd_C=function(knot,g,nu,l,tausq,mk=TRUE){ C.vec=C.eval(knot,g,nu,l,tausq,mk)$vec eval=C.eval(knot,g,nu,l,tausq,mk)$min.eig.val if(eval>0) return(list("cj" = C.vec,"g" = g)) else{ g=g+1 nnd_C(knot,g,nu,l,tausq,mk) } } # computing the eigen values of C using FFT: eigval=function(knot,nu,l,tausq,mk=TRUE){ g=min_g(knot) c.j=nnd_C(knot,g,nu,l,tausq,mk)$cj lambda=Re(fft(c.j)) if(min(lambda)>0) return(lambda) else stop("nnd condition is NOT satisfied!!!") } ################################################################# ########## Samples drawn using Wood and Chan Algorithm ########## ################################################################# samp.WC=function(knot,nu,l,tausq,sseed=1,mk=TRUE){ N=length(knot) lambda=eigval(knot,nu,l,tausq,mk) m=length(lambda) samp.vec=rep(0,N) set.seed(sseed) a=rep(0,m) a[1]=sqrt(lambda[1])*rnorm(1)/sqrt(m) a[(m/2)+1]=sqrt(lambda[(m/2)+1])*rnorm(1)/sqrt(m) i=sqrt(as.complex(-1)) for(j in 2:(m/2)){ uj=rnorm(1); vj=rnorm(1) a[j]=(sqrt(lambda[j])*(uj + i*vj))/(sqrt(2*m)) a[m+2-j]=(sqrt(lambda[j])*(uj - i*vj))/(sqrt(2*m)) } samp=fft(a) samp.vec=Re(samp[1:N]) return(samp.vec) } ############################################# ########## Functions for using ESS ########## ############################################# ESS = function(beta,nu_ess,y,X,sigsq,eta,seeds=1){ thetamin = 0; thetamax = 2*pi; set.seed(seeds) u = runif(1) logy = loglik(y,X,sigsq,eta,beta) + log(u); theta = runif(1,thetamin,thetamax); thetamin = theta - 2*pi; thetamax = theta; betaprime = beta*cos(theta) + nu_ess*sin(theta); while(loglik(y,X,sigsq,eta,betaprime) <= logy){ if(theta < 0) thetamin = theta else thetamax = theta theta = runif(1,thetamin,thetamax) betaprime = beta*cos(theta) + nu_ess*sin(theta) } return(betaprime) } ESS.dec = function(beta,nu_ess,y,X,sigsq,eta,seeds=1){ thetamin = 0; thetamax = 2*pi; set.seed(seeds) u = runif(1) logy = loglik2(y,X,sigsq,eta,beta) + log(u); theta = runif(1,thetamin,thetamax); thetamin = theta - 2*pi; thetamax = theta; betaprime = beta*cos(theta) + nu_ess*sin(theta); while(loglik2(y,X,sigsq,eta,betaprime) <= logy){ if(theta < 0) thetamin = theta else thetamax = theta theta = runif(1,thetamin,thetamax) betaprime = beta*cos(theta) + nu_ess*sin(theta) } return(betaprime) } ## Defining the loglik functions to be used in ESS: ## loglik calculates the log of the likelihood: loglik=function(y,X,sigsq,eta,beta){ mu=y-(X%*%beta) val=eta*sum(beta)-sum(log(1+exp(eta*beta)))-sum(mu^2)/(2*sigsq) return(val) } loglik2=function(y,X,sigsq,eta,beta){ mu=y-(X%*%beta) val=-sum(log(1+exp(eta*beta)))-sum(mu^2)/(2*sigsq) return(val) } ## MH algo for \nu and \ell of Matern kernel: # range.nu=c(0.5,1),range.l=c(0.1,1) : default setting nu.MH2 = function(nu.in,l.in,tau.in,xi.in,knot, range.nu=c(0.5,1),range.l=c(0.1,1), sd.nu=0.1,sd.l=0.1,seed=1){ Kmat = covmat(knot,nu.in,l.in) Linv = inv_chol(Kmat) set.seed(seed) nu.cand = exp(log(nu.in)+rnorm(1,0,sd.nu)) l.cand = exp(log(l.in)+rnorm(1,0,sd.l)) dnu = dunif(nu.cand,range.nu[1],range.nu[2]) dl = dunif(l.cand,range.l[1],range.l[2]) if(dnu > 0 && dl > 0){ Kcand = covmat(knot,nu.cand,l.cand) Linv.cand = inv_chol(Kcand) r = exp(sum(log(diag(Linv.cand)))-sum(log(diag(Linv))))*(nu.cand/nu.in)*(l.cand/l.in) t1 = sum((t(Linv.cand)%*%xi.in)^2) t2 = sum((t(Linv)%*%xi.in)^2) r = r*exp(-(t1 - t2)/(2*tau.in)) alpha = min(r,1) } else{ alpha = 0 Linv.cand = Linv } u = runif(1) nu.out = (u < alpha)*nu.cand + (u >= alpha)*nu.in l.out = (u < alpha)*l.cand + (u >= alpha)*l.in cnt = (u < alpha)*1 + (u >= alpha)*0 L_inv = (u < alpha)*Linv.cand + (u >= alpha)*Linv return(list("nu" = nu.out,"l" = l.out,"count" = cnt,"L_inv"=L_inv)) } ## MH algo for \ell of squared exponential kernel: l.MH = function(l.in,tau.in,xi.in,knot,range.l=c(0.2,1),sd.p=0.05,seed=1){ # corr \in {10^(-6),0.606} for the range of \ell Kmat = covmat.sek(knot,l.in) Linv = inv_chol(Kmat) set.seed(seed) l.cand = exp(log(l.in)+rnorm(1,0,sd.p)) du = dunif(l.cand,range.l[1],range.l[2]) if(du > 0){ Kcand = covmat.sek(knot,l.cand) Linv.cand = inv_chol(Kcand) t1 = sum((t(Linv.cand)%*%xi.in)^2); t2 = sum((t(Linv)%*%xi.in)^2) t3 = sum(log(diag(Linv.cand)))-sum(log(diag(Linv)))-((t1 - t2)/(2*tau.in)) r = (l.cand/l.in)*exp(t3) alpha = min(r,1) }else{ alpha = 0 Linv.cand = Linv } u = runif(1) l.out = (u < alpha)*l.cand + (u >= alpha)*l.in cnt = as.numeric((abs(l.out - l.in) > 0)) L_inv = (u < alpha)*Linv.cand + (u >= alpha)*Linv return(list("l" = l.out,"count" = cnt,"L_inv"=L_inv)) } #################################### #### Error measure calculations #### #################################### err.meas.report=function(data_test,out,reg_i,reg_f,return.mcse){ # data_test: Test data set. # out: The list object that contains the posterior summaries of out-of-sample prediction including # posterior samples of function predictive values; posterior mean and credible intervals of function predictive curves. # reg_i, reg_f: The indicator functions of increasing and flat regions # return.mcse: return.mcse=TRUE returns the effective sample size and standard deviation of the Markov chains # Output: Overall RMSE; RMSE over increasing region; RMSE over the flat region; # Coverage; Effective sample size; Markov chain standard deviation. require(mcmcse) xtest=data_test[,1];yte.tr=data_test[,2]; nte=length(xtest) ftest = yte.tr xte = xtest yte = data_test[,3] fmean=out$fmean;f_low=out$f_low; f_upp=out$f_upp rmse=sqrt(mean((yte.tr-fmean)^2)) if(missing(reg_f)){ rmse.inc="NaN" }else { rmse.flt=sqrt(mean((yte.tr[reg_f(xtest)]-fmean[reg_f(xtest)])^2))} if(missing(reg_i)){ rmse.inc="NaN" }else{ rmse.inc=sqrt(mean((yte.tr[reg_i(xtest)]-fmean[reg_i(xtest)])^2))} cvg=sum((f_low<=yte.tr) & (yte.tr<=f_upp))/nte if(return.mcse){ fhat_sam=out$fhat_sam; ess = mcmcse::ess(t(fhat_sam)) mcse = mcmcse::mcse(fhat_sam)$se }else{ ess="NaN" mcse="NaN" } return(list("rmse"=rmse,"rmse.flt"=rmse.flt,"rmse.inc"=rmse.inc,"cvg"=cvg, "ess"=ess,"mcse"=mcse)) } #### Report the WAIC value of the model on out-of-sample predictions fun.waic <- function(y,fhat_sam,sig_sam){ n = length(y) nmcmc = length(sig_sam) likval = matrix(0,nrow = n, ncol = nmcmc) for(i in 1:n){ for(j in 1:nmcmc){ likval[i,j] = dnorm(y[i], mean = fhat_sam[i,j], sd = sqrt(sig_sam[j])) } } lppd = sum(log(rowMeans(likval))) p_waic1 = 2*sum(log(rowMeans(likval)) - rowMeans(log(likval))) waic1 = -2*(lppd - p_waic1) return("WAIC"=waic1) } #### Report the WAIC value of the model on out-of-sample predictions over a certain region of the function domain fun.waic.reg <- function(x,y,fhat_sam,sig_sam,reg_ind){ yn = y[reg_ind(x)] n = length(yn) nmcmc = length(sig_sam) likval = matrix(0,nrow = n, ncol = nmcmc) for(i in 1:n){ for(j in 1:nmcmc){ likval[i,j] = dnorm(yn[i], mean = fhat_sam[i,j], sd = sqrt(sig_sam[j])) } } lppd = sum(log(rowMeans(likval))) p_waic1 = 2*sum(log(rowMeans(likval)) - rowMeans(log(likval))) waic1 = -2*(lppd - p_waic1) return("WAIC"=waic1) }