R version 2.12.2 (2011-02-25) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-0 Platform: i486-pc-linux-gnu (32-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. > options(digits=6) > # illustrates the use of merge, for merging parameters accross variables: > # Z1=m+e1(s) > # Z2=m+e2(s) > # Z1 and Z2 each have a different variogram, but share the parameter m > # see documentation of gstat() function > library(gstat) Loading required package: sp > d1 = data.frame(x=c(0,2),y=c(0,0),z=c(0,1)) > d2 = data.frame(x=c(0,2),y=c(2,2),z=c(4,5)) > g = gstat(NULL,"d1", z~1,~x+y,d1,model=vgm(1, "Exp", 1)) > g = gstat(g,"d2", z~1,~x+y,d2,model=vgm(1, "Exp", 1), merge=c("d1","d2")) > g = gstat(g, c("d1", "d2"), model = vgm(0.5, "Exp", 1)) > predict.gstat(g, data.frame(x=1,y=1), debug = 32) Intrinsic Correlation found. Good. [using ordinary cokriging] we're at location X: 1 Y: 1 Z: 0 zero block size we're at point X: 1 Y: 1 Z: 0 # X: Matrix: 4 by 1 row 0: 1 row 1: 1 row 2: 1 row 3: 1 [using generalized covariances: max_val - semivariance()] # Covariances (x_i, x_j) matrix C (lower triangle only): Matrix: 4 by 4 row 0: 1 0 0 0 row 1: 0.135335283 1 0 0 row 2: 0.0676676416 0.0295528733 1 0 row 3: 0.0295528733 0.0676676416 0.135335283 1 # X'C-1 X: Matrix: 1 by 1 row 0: 3.24528918 # beta: Vector: dim: 1 2.5 # Cov(beta), (X'C-1 X)-1: Matrix: 1 by 1 row 0: 0.30813895 # Corr(beta): Matrix: 1 by 1 row 0: 1 # X0 (X values at prediction location x0): Matrix: 1 by 2 row 0: 1 1 # BLUE(mu), E(y(x0)) = X0'beta: Vector: dim: 2 2.5 2.5 # Covariances (x_i, x_0), C0: Matrix: 4 by 2 row 0: 0.243116734 0.121558367 row 1: 0.243116734 0.121558367 row 2: 0.121558367 0.243116734 row 3: 0.121558367 0.243116734 # C-1 C0: Matrix: 4 by 2 row 0: 0.206482174 0.089386867 row 1: 0.206482174 0.089386867 row 2: 0.089386867 0.206482174 row 3: 0.089386867 0.206482174 # [a] Cov_ij(B,B) or Cov_ij(0,0): Matrix: 2 by 2 row 0: 1 0.5 row 1: 0.5 1 # [c] (x0-X'C-1 c0)'(X'C-1 X)-1(x0-X'C-1 c0): Matrix: 2 by 2 row 0: 0.0513599204 0.0513599204 row 1: 0.0513599204 0.0513599204 # [b] c0'C-1 c0: Matrix: 2 by 2 row 0: 0.122129987 0.0936621582 row 1: 0.0936621582 0.122129987 # Best Linear Unbiased Predictor: Vector: dim: 2 2.03161877 2.96838123 # MSPE ([a]-[b]+[c]): Matrix: 2 by 2 row 0: 0.929229934 0.457697762 row 1: 0.457697762 0.929229934 # kriging weights: Matrix: 4 by 2 row 0: 0.308547653 0.191452347 row 1: 0.308547653 0.191452347 row 2: 0.191452347 0.308547653 row 3: 0.191452347 0.308547653 x y d1.pred d1.var d2.pred d2.var cov.d1.d2 1 1 1 2.03162 0.92923 2.96838 0.92923 0.457698 > > # Z1 and Z2 share a regression slope: > g = gstat(NULL,"d1", z~x,~x+y,d1,model=vgm(1, "Exp", 1)) > g = gstat(g,"d2", z~x,~x+y,d2,model=vgm(1, "Exp", 1), + merge=list(c("d1",2,"d2",2))) > g = gstat(g, c("d1", "d2"), model = vgm(0.5, "Exp", 1)) > predict.gstat(g, data.frame(x=1,y=1), debug = 32) Intrinsic Correlation found. Good. [using universal cokriging] we're at location X: 1 Y: 1 Z: 0 zero block size we're at point X: 1 Y: 1 Z: 0 # X: Matrix: 4 by 3 row 0: 1 0 0 row 1: 1 2 0 row 2: 0 0 1 row 3: 0 2 1 [using generalized covariances: max_val - semivariance()] # Covariances (x_i, x_j) matrix C (lower triangle only): Matrix: 4 by 4 row 0: 1 0 0 0 row 1: 0.135335283 1 0 0 row 2: 0.0676676416 0.0295528733 1 0 row 3: 0.0295528733 0.0676676416 0.135335283 1 # X'C-1 X: Matrix: 3 by 3 row 0: 1.77460693 1.62264459 -0.151962334 row 1: 1.62264459 7.67605004 1.62264459 row 2: -0.151962334 1.62264459 1.77460693 # beta: Vector: dim: 3 2.10985743e-16 0.5 4 # Cov(beta), (X'C-1 X)-1: Matrix: 3 by 3 row 0: 0.793362513 -0.225694871 0.274305129 row 1: -0.225694871 0.225694871 -0.225694871 row 2: 0.274305129 -0.225694871 0.793362513 # Corr(beta): Matrix: 3 by 3 row 0: 1 -0.533365606 0.34575005 row 1: -0.533365606 1 -0.533365606 row 2: 0.34575005 -0.533365606 1 # X0 (X values at prediction location x0): Matrix: 3 by 2 row 0: 1 0 row 1: 1 1 row 2: 0 1 # BLUE(mu), E(y(x0)) = X0'beta: Vector: dim: 2 0.5 4.5 # Covariances (x_i, x_0), C0: Matrix: 4 by 2 row 0: 0.243116734 0.121558367 row 1: 0.243116734 0.121558367 row 2: 0.121558367 0.243116734 row 3: 0.121558367 0.243116734 # C-1 C0: Matrix: 4 by 2 row 0: 0.206482174 0.089386867 row 1: 0.206482174 0.089386867 row 2: 0.089386867 0.206482174 row 3: 0.089386867 0.206482174 # [a] Cov_ij(B,B) or Cov_ij(0,0): Matrix: 2 by 2 row 0: 1 0.5 row 1: 0.5 1 # [c] (x0-X'C-1 c0)'(X'C-1 X)-1(x0-X'C-1 c0): Matrix: 2 by 2 row 0: 0.20356416 -0.100844319 row 1: -0.100844319 0.20356416 # [b] c0'C-1 c0: Matrix: 2 by 2 row 0: 0.122129987 0.0936621582 row 1: 0.0936621582 0.122129987 # Best Linear Unbiased Predictor: Vector: dim: 2 0.5 4.5 # MSPE ([a]-[b]+[c]): Matrix: 2 by 2 row 0: 1.08143417 0.305493523 row 1: 0.305493523 1.08143417 # kriging weights: Matrix: 4 by 2 row 0: 0.5 6.9388939e-17 row 1: 0.5 -9.71445147e-17 row 2: -4.16333634e-17 0.5 row 3: 2.77555756e-17 0.5 x y d1.pred d1.var d2.pred d2.var cov.d1.d2 1 1 1 0.5 1.08143 4.5 1.08143 0.305494 >