[R] Simulating scalar-valued stationary Gaussian processes
Ben Bolker
bolker at ufl.edu
Sat May 13 14:07:48 CEST 2006
Mikael Anderson <mikael.anderson <at> gmail.com> writes:
>
> Hi,
>
> I have a sample of size 100 from a function in interval [0,1] which can be
> assumed to come from a scalar-valued stationary Gaussian process. There are
> about 500 observation points in the interval. I need an effective and fast
> way to simulate from the Gaussian process conditioned on the available data.
> I can of course estimate the mean and 500x500 covariance matrix from data.
> I have searched both the Rsite and Internet in general but have not found
> any function or program which could work directly for this type of data.
> There are some packages in R to simulate from Gaussian processes(e.g. geoR,
> waveslim and spectralGP) but I couldn't figure out if any function in these
> packages could work in my framework.
>
> Given the size of covariance matrix I was thinking about using SVD to reduce
> the size but I would like to hear from the list members if there are more
> effective ways of doing this.
>
One option would be to Fourier transform, randomize the
phase, and back-transform: something like
function (x, covar = NULL, argtype = "even", intype = "cov")
{
l <- length(x)
if (!is.null(covar))
x <- x * sqrt(covar/var(x))
y <- fft(x)
if (argtype == "even") {
if (l%%2 != 0) {
cat("length must be even, for now\n")
return(NA)
}
arg0 <- runif(l/2 - 1, -pi, pi)
arg1 <- c(0, arg0, 0, -rev(arg0))
}
else if (argtype == "zero")
arg1 <- rep(0, length(y))
else if (argtype == "random")
arg1 <- runif(l, -pi, pi)
z <- complex(modulus = Mod(y), argument = arg1)
y <- fft(z, inverse = TRUE)/l
chop(y)
}
of course this also depends exactly what you mean by
"conditioned on the available data" -- and how do you
feel about parametric models? RandomFields has lots of
tools for conditional simulation, if you don't mind fitting
a parametric covariance model to your data ...
Ben Bolker
More information about the R-help
mailing list