R/LKrig.sim.conditional.R

In LatticeKrig: Multi-Resolution Kriging Based on Markov Random Fields

# LatticeKrig  is a package for analysis of spatial data written for
# the R software environment .
# Copyright (C) 2016
# University Corporation for Atmospheric Research (UCAR)
# Contact: Douglas Nychka, nychka@ucar.edu,
# National Center for Atmospheric Research, PO Box 3000, Boulder, CO 80307-3000
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with the R software environment if not, write to the Free Software
# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
# or see http://www.r-project.org/Licenses/GPL-2

LKrig.sim.conditional <- function(LKrigObj, M = 1, x.grid = NULL, 
    grid.list = NULL, nx = 80, ny = 80, ..., Z.grid = NULL, seed=42, verbose=FALSE) {
    # generate grid if not specified
    if (is.null(x.grid)) {
        if (is.null(grid.list)) {
            grid.list <- fields.x.to.grid(LKrigObj$x, nx = nx, ny = ny)
        }
        x.grid <- make.surface.grid(grid.list)
    }
    # coerce x.grid to a matrix.  
    x.grid<- as.matrix( x.grid)
    if( verbose){
    	cat("LKrig.sim.conditional: x.grid")
    	print( x.grid)
    }
    # NOTE: the name x.grid may be misleading because it just needs to a  matrix with rows
    # indexing th locations and need not follow any regular pattern.
    # Now generate the error surfaces
    # begin block
    # create vector of seeds if needed
    if( length(seed)==1){
        seeds<- seed + ( 0:(M-1))
      }
    #
    g.conditional.draw <-    matrix(NA, ncol = M, nrow = nrow(x.grid))

    N <- nrow(LKrigObj$y)  
    # complete set of locations to evaluate the field must
    # include the observations too
    PHIGrid<- LKrig.basis(x.grid,LKrigObj$LKinfo) 
    if( verbose){
    	cat("LKrig.sim.conditional: dim(PHIGrid)",  dim(PHIGrid), fill=TRUE)
    }
   
    # predicted field at grid from the actual data
    spatialPart<- (PHIGrid%*% LKrigObj$c.coef)
    ghat <-  spatialPart
    if( !is.null(LKrigObj$LKinfo$fixedFunction) ){
      d.coef.draw<- matrix(NA, ncol = M, nrow = length( LKrigObj$d.coef) )
   }
   else{
   	d.coef.draw<- matrix(NA, ncol = M, nrow=1)
   	}
   	ghat<- predict( LKrigObj, x= x.grid, Z = Z.grid )
    for (k in 1:M) {
        cat(k, " ") # numbers to entertain ...
        out<- simConditionalDraw( k, LKrigObj, ghat, x.grid, Z.grid,  PHIGrid,
                                 seeds, ..., verbose=verbose)
        if( !is.null(LKrigObj$LKinfo$fixedFunction) ){
        d.coef.draw[,k] <- out$d.coef
        }
        g.conditional.draw[, k] <- out$g.conditional
    }
   	cat(k, " ", fill=TRUE)
    # sample standard deviation across ensemble members
   	SE<-  apply(g.conditional.draw, 1, FUN=sd)
    return(list(x.grid = x.grid, ghat = ghat, g.draw = g.conditional.draw,
                           SE = SE,
                           d.coef.draw= d.coef.draw))
}

simConditionalDraw <- function(index=1,  LKrigObj, ghat, x.grid, Z.grid,
                               PHIGrid, seeds= 123,  verbose=FALSE){
require(LatticeKrig)
        set.seed( seeds[index] )
# generate process at grid and also on the observation locations.
        simCoefficients<- LKrig.sim(LKinfo = LKrigObj$LKinfo, just.coefficients=TRUE)
        g.unconditional.data <-LKrigObj$wX %*%simCoefficients 
        g.unconditional.data <- sqrt(LKrigObj$rho.MLE) * g.unconditional.data
        g.unconditional.grid <-sqrt(LKrigObj$rho.MLE) *PHIGrid%*%simCoefficients 
        # generate a synthetic data set with fixed part set to zero.
        N<- length( LKrigObj$y)
        y.synthetic.data <- g.unconditional.data + LKrigObj$sigma.MLE * 
            rnorm(N)
# this may confusing. divide by  sqrt(weights) to cancel out this term in the
# wX matrix for  g.unconditional.data and to adjust measurement error variance         
        y.synthetic.data<- y.synthetic.data / sqrt(LKrigObj$weights)
        # use LKrig to find the predictions for the xgrid locations
        # NOTE that LKrig will still estimate the fixed part.
        # and it is important to include this part of estimate
        if(verbose){
           cat("simConditionalDraw Call to LKrig: ", fill=TRUE)
        }
        obj.fit.synthetic <- LKrig(LKrigObj$x, y.synthetic.data,
                                   LKinfo = LKrigObj$LKinfo,
                                       wX = LKrigObj$wX,
                                       wU = LKrigObj$wU,
                                   lambda = LKrigObj$lambda,
                                        Z = LKrigObj$Z,
                                  weights = LKrigObj$weights,
                             use.cholesky = LKrigObj$Mc, 
                                  verbose = verbose)      
        #
        # predict field
        spatialPart<- (PHIGrid%*% obj.fit.synthetic$c.coef)
        ghat.synthetic<-  spatialPart
        if( !is.null(LKrigObj$LKinfo$fixedFunction) ){       	
                 fixedPart<- predict(
                     obj.fit.synthetic, xnew=x.grid, Znew = Z.grid,
                     just.fixed=TRUE)
                 d.coef <- obj.fit.synthetic$d.coef
                 ghat.synthetic<- ghat.synthetic + fixedPart
         }
         else{
         	d.coef<- NA
         	}
         	
        # add prediction error to the condition mean from the actual data
        g.conditional <- ghat + (g.unconditional.grid -  ghat.synthetic)
        # NOTE: sampling variablity for fixed part is built in too
        # because d.coef are estimated and included in the prediction. 
return(
       list( g.conditional = g.conditional, d.coef = d.coef) )
}