R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Functions for simulating a multiresolution process following...

LKrig.sim

R Documentation

Functions for simulating a multiresolution process following the Lattice Krig covariance model.

Description

The fields are Gaussian and can be either simulated unconditionally or conditional on the field values and a set of irregular locations.

Usage

#
LKrig.sim(x1, LKinfo, M=1,just.coefficients = FALSE)
LKrig.sim.conditional( LKrigObj, M=1, x.grid= NULL, grid.list=NULL,
                           nx=80, ny=80,...,Z.grid=NULL, seed=42, verbose=FALSE)

simConditionalDraw(index = 1, LKrigObj, ghat, x.grid, Z.grid, PHIGrid,
                 seeds = 123, verbose = FALSE)

Arguments

`grid.list`	Specifies a grid of spatial locations using the `grid.list` format ( help(grid.list)). These are the locations used to evaluate the fields generated from conditional simulation. The default is to generate an 80X80 grid based on range of the observations.
`just.coefficients`	If TRUE just simulates the coefficients from the Markov Random field.
`LKinfo`	A list with components that give the information describing a multiresolution basis with a Markov random field used for the covariance of the basis coefficients. This list is created in `LKrig` and is returned as part of the output object or in a more hands on manner directly using `LKrigSetup` (See section on returned value below for this list's description.)
`M`	Number of independent simulated fields.
`nx`	Number of grid points in x coordinate for output grid.
`ny`	Number of grid points in y coordinate for output grid.
`LKrigObj`	An `LKrig` object, i.e. the output list returned by `LKrig`.
`seed`	Seed to set random number generator.
`x1`	A two column matrix of 2-dimension locations to evaluate basis functions or the first set of locations to evaluate the covariance function or the locations for the simulated process. Rows index the different locations: to be precise `x1[i,1:2]` are the "x" and "y" coordinates for the `i` th location.
`x.grid`	Locations to evaluate conditional fields. This is in the form of a two column matrix where each row is a spatial location.
`Z.grid`	The covariates that are associated with the x.grid values. This is useful for conditional simulation where the fields are evaluated at x.grid locations and using covariate values Z.grid. Z.grid is matrix with columns indexing the different covariates and rows indexed by the x.grid locations.
`...`	Arguments to be passed to the LKrig function to specify the spatial estimate. These are components in addition to what is in the LKinfo list returned by LKrig.
`verbose`	If TRUE prints out debugging information.
`ghat`	The predicted surface at the grid.
`index`	The index for the random seed to use in the vector `seeds`.
`PHIGrid`	Basis function matrix at grid points.
`seeds`	A vector of random seeds.

Details

The simulation of the unconditional random field is done by generating a draw from the multi-resolution coefficients using a Cholesky decomposition and then multiplying by the basis functions to evaluate the field at arbitrary points. Currently, there is no provision to exploit the case when one wants to simulate the field on a regular grid. The conditional distribution is a draw from the multivariate normal for the random fields conditioned on the observations and also conditioned on covariance model and covariance parameters. If the nugget/measurement error variance is zero then any draw from the conditional distribution will be equal to the observations at the observation locations. In the past conditional simulation was known to be notoriously compute intensive, but the major numerical problems are finessed here by exploiting sparsity of the coefficient precision matrix.

The conditional field is found using a simple trick based on the linear statistics for the multivariate normal. One generates an unconditional field that includes the field values at the observations. From this realization one forms a synthetic data set and uses LKrig to predict the remaining field based on the synthetic observations. The difference between the predicted field and the realization (i.e. the true field) is a draw from the conditional distribution with the right covariance matrix. Adding the conditional mean to this result one obtains a draw from the full conditional distribution. This algorithm can also be interpreted as a variant on the bootstrap to determine the estimator uncertainty. The fixed part of the model is also handled correctly in this algorithm. See the commented source for LKrig.sim.conditional for the details of this algorithm.

simConditionalDraw is low level function that is called to generate each ensemble member i.e. each draw from the conditional distribution. The large number of arguments is to avoid recomputing many common elements during the loop in generating these draws. In particular passing the basis function matrices avoid having to recompute the normalization at each step, often an intensive computation for a large grid.

Value

LKrig.sim: A matrix with dimensions of nrow(x1) by M of simulated values at the locations x1.

LKrig.sim.conditional: A list with the components.

xgrid: The locations where the simulated field(s) are evaluated.
ghat: The conditional mean at the xgrid locations.
g.draw: A matrix with dimensions of nrow(x.grid) by M with each column being an independent draw from the conditional distribution.

Author(s)

Doug Nychka

Examples

# Load ozone data set
  data(ozone2)  
  x<-ozone2$lon.lat
  y<- ozone2$y[16,]
# Find location that are not 'NA'.
# (LKrig is not set up to handle missing observations.)
  good <-  !is.na( y)
  x<- x[good,]
  y<- y[good]
  LKinfo<- LKrigSetup( x,NC=20,nlevel=1, alpha=1, lambda= .3 , a.wght=5)
# BTW lambda is close to MLE 
# Simulating this  LKrig process
# simulate 4 realizations of process and plot them
# (these have unit marginal variance)
  xg<- make.surface.grid(list( x=seq( -87,-83,,40), y=seq(36.5, 44.5,,40)))
  out<- LKrig.sim(xg, LKinfo,M=4)
## Not run: 
  set.panel(2,2)
  for( k in 1:4){
    image.plot( as.surface( xg, out[,k]), axes=FALSE) }

## End(Not run)
  obj<- LKrig(x,y,LKinfo=LKinfo)
  O3.cond.sim<- LKrig.sim.conditional( obj, M=3,nx=40,ny=40) 
## Not run: 
  set.panel( 2,2)
  zr<- range( c(  O3.cond.sim$draw,  O3.cond.sim$ghat), na.rm=TRUE)
  coltab<- tim.colors()
  image.plot( as.surface( O3.cond.sim$x.grid, O3.cond.sim$ghat), zlim=zr)
  title("Conditional mean")
  US( add=TRUE)
  for( k in 1:3){
    image( as.surface( O3.cond.sim$x.grid, O3.cond.sim$g.draw[,k]),
              zlim=zr, col=coltab)
    points( obj$x, cex=.5)
    US( add=TRUE)
  }
  set.panel()

## End(Not run)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(LatticeKrig)
Loading required package: spam
Loading required package: grid
Spam version 1.3-0 (2015-10-24) is loaded.
Type 'help( Spam)' or 'demo( spam)' for a short introduction 
and overview of this package.
Help for individual functions is also obtained by adding the
suffix '.spam' to the function name, e.g. 'help( chol.spam)'.

Attaching package: 'spam'

The following objects are masked from 'package:base':

    backsolve, forwardsolve

Loading required package: fields
Loading required package: maps

 # maps v3.1: updated 'world': all lakes moved to separate new #
 # 'lakes' database. Type '?world' or 'news(package="maps")'.  #


> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/LatticeKrig/LKrig.sim.Rd_%03d_medium.png", width=480, height=480)
> ### Name: LKrig.sim
> ### Title: Functions for simulating a multiresolution process following the
> ###   Lattice Krig covariance model.
> ### Aliases: LKrig.sim LKrig.sim.conditional simConditionalDraw
> ### Keywords: spatial
> 
> ### ** Examples
> 
> # Load ozone data set
>   data(ozone2)  
>   x<-ozone2$lon.lat
>   y<- ozone2$y[16,]
> # Find location that are not 'NA'.
> # (LKrig is not set up to handle missing observations.)
>   good <-  !is.na( y)
>   x<- x[good,]
>   y<- y[good]
>   LKinfo<- LKrigSetup( x,NC=20,nlevel=1, alpha=1, lambda= .3 , a.wght=5)
> # BTW lambda is close to MLE 
> # Simulating this  LKrig process
> # simulate 4 realizations of process and plot them
> # (these have unit marginal variance)
>   xg<- make.surface.grid(list( x=seq( -87,-83,,40), y=seq(36.5, 44.5,,40)))
>   out<- LKrig.sim(xg, LKinfo,M=4)
> ## Not run: 
> ##D   set.panel(2,2)
> ##D   for( k in 1:4){
> ##D     image.plot( as.surface( xg, out[,k]), axes=FALSE) }
> ## End(Not run)
>   obj<- LKrig(x,y,LKinfo=LKinfo)
>   O3.cond.sim<- LKrig.sim.conditional( obj, M=3,nx=40,ny=40) 
1  2  3  > ## Not run: 
> ##D   set.panel( 2,2)
> ##D   zr<- range( c(  O3.cond.sim$draw,  O3.cond.sim$ghat), na.rm=TRUE)
> ##D   coltab<- tim.colors()
> ##D   image.plot( as.surface( O3.cond.sim$x.grid, O3.cond.sim$ghat), zlim=zr)
> ##D   title("Conditional mean")
> ##D   US( add=TRUE)
> ##D   for( k in 1:3){
> ##D     image( as.surface( O3.cond.sim$x.grid, O3.cond.sim$g.draw[,k]),
> ##D               zlim=zr, col=coltab)
> ##D     points( obj$x, cex=.5)
> ##D     US( add=TRUE)
> ##D   }
> ##D   set.panel()
> ## End(Not run)
> 
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>