1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
#
# fields is a package for analysis of spatial data written for
# the R software environment.
# Copyright (C) 2024 Colorado School of Mines
# 1500 Illinois St., Golden, CO 80401
# Contact: Douglas Nychka, douglasnychka@gmail.com,
#
# 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
##END HEADER
sim.spatialProcess<- function(object, xp, M = 1,
verbose = FALSE, ...) {
# important variance parameters estimated from the data
tau2 <- (object$summary["tau"])^2
sigma2 <- (object$summary["sigma2"])
xp<- as.matrix( xp)
#
# check for unique rows of data locations
if( any(duplicated(object$x)) ){
stop("data locations should be unique")
}
#
# find sizes of arrays
m <- nrow(xp)
n<- nrow( object$x)
N <- length(object$y)
if (verbose) {
cat("m,n,N", m,n, N, fill = TRUE)
}
#
# find indices of all rows of xp that correspond to rows of
# xM and then collapse x to unique rows.
if( any( duplicated(object$x)) ){
stop('Can not handle replicated locations')}
if( any( duplicated(xp)) ){
stop('Can not handle repeated
prediction locations')}
#
x<- as.matrix(rbind( object$x, xp))
rep.x.info <- fields.duplicated.matrix(x)
# find uniuqe locations.
ind<- !duplicated(rep.x.info)
xUnique <- as.matrix(x[ind, ])
if (verbose) {
cat('full x and predicted locations without duplicates ', fill = TRUE)
print(xUnique)
}
N.full <- nrow(xUnique)
if (verbose) {
cat("N.full", N.full, fill = TRUE)
}
if( N.full > 5000){
cat("WARNING: Number of locations for conditional simulation is large ( >5000)
this may take some time to compute or exhaust the memory.",
fill=FALSE)
}
# these give locations in x matrix to reconstruct xp matrix
xp.ind <- rep.x.info[(1:m) + n]
if (verbose) {
print(N.full)
print(xUnique)
}
if (verbose) {
cat("reconstruction of xp from collapsed locations",
fill = TRUE)
print(xUnique[xp.ind, ])
}
#
# Sigma is full covariance at the data locations and at prediction points.
# not to be confused with the lowercase tau that is the nugget variance
#
Sigma <- sigma2 * do.call(object$cov.function.name, c(object$args,
list(x1 = xUnique, x2 = xUnique)))
#
# square root of Sigma for simulating field
# Cholesky is fast but not very stable.
#
# the following code line is similar to chol(Sigma)-> Scol
# but adds possible additional arguments controlling the Cholesky
# from the mKrig object.
# x has has been winnowed down to unique rows so that
# Sigma has full rank.
#
Schol <- do.call("chol", c(list(x = Sigma), object$chol.args))
#
# output matrix to hold results
out <- matrix(NA, ncol = M, nrow = m)
#
# find conditional mean field from initial fit
# (these are added at the predict step).
#
h.hat <- predict(object, xnew=xp, ...)
#
# NOTE: fixed part of model (null space) need not be simulated
# because the estimator is unbiased for this part.
# the variability is still captured because the fixed part
# is still estimated as part of the predict step below
# create synthetic data
for (k in 1:M) {
if( k%%10==1){
cat( k," ")
}
# simulate full field
h <- t(Schol) %*% rnorm(N.full)
# value of simulated field at observations
h.data <- h[1:n]
#
y.synthetic <- h.data + sqrt(tau2/object$weights)*rnorm(n)
# predict at xp using these data
# and subtract from 'true' value,
# note that true values of field have to be expanded in the
# case of common locations between object$x and xp.
h.true <- (h[xp.ind])
# temp.error <- predict(object, xnew=xp, ynew = y.synthetic,
# Z=Zp, ...) - h.true
temp.error <- predict(object, xnew=xp, ynew = y.synthetic,
...) - h.true
# add the error to the actual estimate (conditional mean)
out[,k ] <- h.hat + temp.error
}
out
}
|
