R/metropolis.krige.R

In lbaole17/krige.ext: Geospatial Kriging with Metropolis Sampling (Experimental Version)

#SAMPLER
metropolis.krige <- function(y,X,east,north,powered.exp=2,mcmc.samples=100,spatial.share=0.5,range.share=0.5,beta.var=10,
	range.tol=0.05,b.tune=1.0,nugget.tune=10.0,psill.tune=1.0,
	message = NULL,save = NULL,time = NULL,checks = TRUE){

  #ERROR CHECKS
	if(checks==TRUE){
	if(any(is.na(y))) stop("missing values in 'y' ")
	if(TRUE %in% apply(X, 2, function(x) any(is.na(x)))) stop("missing values in 'X' ")
	if(any(is.na(east))) stop("missing values in 'east' ")
	if(any(is.na(north))) stop("missing values in 'north' ")
  if(powered.exp<=0 | powered.exp>2) stop("powered.exp must be greater than 0 and less than or equal to 2.")
  if(spatial.share<=0 | spatial.share>=1) stop("spatial.share must be between 0 and 1.")
  if(range.share<=0) stop("range.share must be greater than 0.")
  if(range.tol<=0 | range.tol>=1) stop("p.range.tol must be between 0 and 1.")
  if(beta.var<=0) stop("beta.var must be greater than 0.")
  if(b.tune<=0) stop("b.tune must be greater than 0.")
  if(nugget.tune<=0) stop("nugget.tune must be greater than 0.")
  if(psill.tune<=0) stop("psill.tune must be greater than 0.")
}

  #INTERNAL FUNCTION, SIMPLIFIED VERSION OF mvrnorm FROM MASS
  simp.mvrnorm<-function(n=1, mu, Sigma){
    p <- length(mu)
    eS <- eigen(Sigma, symmetric = TRUE)
    ev <- eS$values
    X <- matrix(rnorm(p * n), n)
    X <- drop(mu)+eS$vectors%*%diag(sqrt(pmax(ev,0)),p)%*%t(X)
    if (n == 1)
        drop(X)
    else t(X)
    }

  #DEFINE STANDING PARAMETERS
  dist.mat<-k_distmat(cbind(east,north))
  max.distance<-max(dist.mat)
  min.distance<-min(dist.mat[dist.mat>0])
  init.ols<-lm(y~X-1)
  b.cand.var<-b.tune*vcov(init.ols)
  err.var<-var(resid(init.ols))
  beta.A<-round(100*range.share/((-log(range.tol))^(1/powered.exp)))
  beta.B<-round(100*(1-(range.share/((-log(range.tol))^(1/powered.exp)))))

  # STARTING VALUES
  mcmc.mat<-matrix(NA,nrow=mcmc.samples,ncol=3+ncol(X))
  if(is.null(dimnames(X)[[2]])) {
  	mynames<-rep(".",dim(X)[2])
  	for(k in 1:dim(X)[2]) mynames[k]<-paste("V",k,sep="")
  	dimnames(X)[[2]]<-mynames
  	}
  dimnames(mcmc.mat)[[2]]<-c("tau2","phi","sigma2",dimnames(X)[[2]])
  n.vars<-ncol(X)
  accepted.beta<-0
  accepted.nugget<-0
  accepted.decay<-0
  accepted.psill<-0
  start.decay<-((-log(range.tol))^(1/powered.exp))/(max(dist.mat)*range.share)
  mcmc.mat[1,]<-c((1-spatial.share)*err.var,start.decay,spatial.share*err.var,init.ols$coef)
  local.Sigma<-NULL

	time0 <- Sys.time()
  #START OF MARKOV CHAIN
  for (i in 1:(nrow(mcmc.mat)-1))  {
		# Message functions
		if (!is.null(save)) if (i%%save == 0) {
			save(mcmc.mat, file = "mcmc-mat.Rdata")
			cat(i, " iterations saved", ".\n", sep = "")}
		if (!is.null(time)){
			if (i %in% time){
				time1 <- Sys.time()
				time.diff <- time1 - time0
				est.time <- (time.diff*(mcmc.samples-i))/i
				if (as.numeric(time.diff, units = "mins") > 120) {t.unit1 <- "hours"} else {t.unit1 <- "mins" }
				if (as.numeric(est.time, units = "mins") > 120) {t.unit2 <- "hours"} else {t.unit2 <- "mins" }
				cat(i, " iterations ", "complete in ", round(as.numeric(time.diff, units = t.unit1), 2), " ", t.unit1, "; ",
				"the rest will take roughly ", round(as.numeric(est.time, units = t.unit2), 2), " ", t.unit2,".\n", sep = "")}
			}
		if (is.null(message)) message <- mcmc.samples/10
    if (i%%message == 0) {
      cat("Iteration ", i, ": Nugget=", mcmc.mat[i, 1], ". Decay=", mcmc.mat[i, 2],
          ". Partial sill=", mcmc.mat[i, 3], ".\n", sep = "")}
		if (i%%message == 0) {
			beta.rate.temp <-round(100*accepted.beta/i)
			tau2.rate.temp <-round(100*accepted.nugget/i)
			phi.rate.temp <-round(100*accepted.decay/i)
			sigma2.rate.temp <-round(100*accepted.psill/i)
			cat("Iteration ", i, ": Acceptance percentages: Coefficients=",beta.rate.temp,"%. Nugget=",tau2.rate.temp,
						"%. Decay=",phi.rate.temp,"%. Partial sill=",sigma2.rate.temp,"%.", ".\n", sep="")}

    old.beta<-mcmc.mat[i,4:ncol(mcmc.mat)]
    old.var<-var(y-X%*%old.beta)
    beta<-simp.mvrnorm(n=1,mu=old.beta,Sigma=b.cand.var)
    local.share<-mcmc.mat[i,3]/(mcmc.mat[i,1]+mcmc.mat[i,3])


    local.tau2.shape<-1+nugget.tune/(1-local.share)
    local.sigma2.shape<-1+psill.tune/local.share
    tau2<-1/rgamma(1,shape=local.tau2.shape,rate=nugget.tune*old.var)
    sigma2<-1/rgamma(1,shape=local.sigma2.shape,rate=psill.tune*old.var)


    phi<-1/(max.distance*rbeta(1,shape1=beta.A,shape2=beta.B))
    if(is.null(local.Sigma)) local.Sigma<-ifelse(dist.mat>0, mcmc.mat[i,3]*exp(-abs(mcmc.mat[i,2]*dist.mat)^powered.exp), mcmc.mat[i,1]+mcmc.mat[i,3])
    current <- krige.posterior(mcmc.mat[i,1],mcmc.mat[i,2],mcmc.mat[i,3],old.beta,
    		y,X,east,north,semivar.exp=powered.exp,p.spatial.share=spatial.share,
    			p.range.share=range.share,p.range.tol=range.tol,p.beta.var=beta.var,tot.var=err.var,
					local.Sigma=local.Sigma,distance = dist.mat,max.distance=max.distance,checks=checks)
    candidate.beta <- krige.posterior(mcmc.mat[i,1],mcmc.mat[i,2],mcmc.mat[i,3],beta,
    		y,X,east,north,semivar.exp=powered.exp,p.spatial.share=spatial.share,
    			p.range.share=range.share,p.range.tol=range.tol,p.beta.var=beta.var,tot.var=err.var,
					local.Sigma=local.Sigma,distance = dist.mat,max.distance=max.distance,checks=checks)
    candidate.nugget <- krige.posterior(tau2,mcmc.mat[i,2],mcmc.mat[i,3],old.beta,
    		y,X,east,north,semivar.exp=powered.exp,p.spatial.share=spatial.share,
    			p.range.share=range.share,p.range.tol=range.tol,p.beta.var=beta.var,tot.var=err.var,
					local.Sigma=NULL,distance = dist.mat,max.distance=max.distance,checks=checks)
    candidate.decay <- krige.posterior(mcmc.mat[i,1],phi,mcmc.mat[i,3],old.beta,
    		y,X,east,north,semivar.exp=powered.exp,p.spatial.share=spatial.share,
    			p.range.share=range.share,p.range.tol=range.tol,p.beta.var=beta.var,tot.var=err.var,
					local.Sigma=NULL,distance = dist.mat,max.distance=max.distance,checks=checks)
    candidate.psill <- krige.posterior(mcmc.mat[i,1],mcmc.mat[i,2],sigma2,old.beta,
    		y,X,east,north,semivar.exp=powered.exp,p.spatial.share=spatial.share,
    			p.range.share=range.share,p.range.tol=range.tol,p.beta.var=beta.var,tot.var=err.var,
					local.Sigma=NULL,distance = dist.mat,max.distance=max.distance,checks=checks)
	a.beta<-exp(candidate.beta-current)
	a.nugget<-exp(candidate.nugget-current)
    a.decay<-exp(candidate.decay-current)
	a.psill<-exp(candidate.psill-current)
    #if (is.na(a.beta) || is.na(a.nugget) || is.na(a.decay) || is.na(a.psill)) {print(paste("Undefined acceptance ratio at iteration ",i)); a.beta<-a.nugget<-a.decay<-a.psill<-0}
    if (a.beta > runif(1)) {
      accepted.beta <- accepted.beta + 1
      mcmc.mat[(i+1),4:ncol(mcmc.mat)] <- beta
    	  } else mcmc.mat[(i+1),4:ncol(mcmc.mat)] <- mcmc.mat[i,4:ncol(mcmc.mat)]
    if (a.nugget > runif(1)) {
      accepted.nugget<-accepted.nugget + 1
      mcmc.mat[(i+1),1] <- tau2
      local.Sigma<-NULL
      } else mcmc.mat[(i+1),1] <- mcmc.mat[i,1]
    if (a.decay > runif(1)) {
      accepted.decay<-accepted.decay+1
      mcmc.mat[(i+1),2] <- phi
      local.Sigma<-NULL
      } else mcmc.mat[(i+1),2] <- mcmc.mat[i,2]
    if (a.psill > runif(1)) {
      accepted.psill<-accepted.psill + 1
      mcmc.mat[(i+1),3] <- sigma2
      local.Sigma<-NULL
      } else mcmc.mat[(i+1),3] <- mcmc.mat[i,3]
  }
  beta.rate<-round(100*accepted.beta/(nrow(mcmc.mat)-1))
  tau2.rate<-round(100*accepted.nugget/(nrow(mcmc.mat)-1))
  phi.rate<-round(100*accepted.decay/(nrow(mcmc.mat)-1))
  sigma2.rate<-round(100*accepted.psill/(nrow(mcmc.mat)-1))
  print(paste("Acceptance percentages: Coefficients=",beta.rate,"%. Nugget=",tau2.rate,"%. Decay=",phi.rate,"%. Partial sill=",sigma2.rate,"%.",sep=""))
  if(beta.rate<20) print("Coefficient acceptance rate is low, which can indicate slow mixing. You may need to run many iterations or consider adjusting the tuning parameter *b.tune* to increase acceptance.")
  if(beta.rate>70) print("Coefficient acceptance rate is high. If coefficients are nonconvergent, consider adjusting the tuning parameter *b.tune* to reduce acceptance.")
  if(tau2.rate<20) print("Nugget acceptance rate is low, which can indicate slow mixing. You may need to run many iterations or consider adjusting the tuning parameter *nugget.tune* to increase acceptance.")
  if(tau2.rate>70) print("Nugget acceptance rate is high. If tau2 is nonconvergent, consider adjusting the tuning parameter *nugget.tune* to reduce acceptance.")
  if(phi.rate<20) print("Decay term acceptance rate is low, which can indicate slow mixing. You may need to run many iterations or consider adjusting the tolerance parameter *range.tol* to increase acceptance.")
  if(phi.rate>70) print("Decay term acceptance rate is high. If phi is nonconvergent, consider adjusting the tolerance parameter *range.tol* to reduce acceptance.")
  if(sigma2.rate<20) print("Partial sill acceptance rate is low, which can indicate slow mixing. You may need to run many iterations or consider adjusting the tuning parameter *psill.tune* to increase acceptance.")
  if(sigma2.rate>70) print("Partial sill acceptance rate is high. If sigma2 is nonconvergent, consider adjusting the tuning parameter *psill.tune* to reduce acceptance.")
  return(mcmc.mat)
}