[R] Simulating scalar-valued stationary Gaussian processes

[Ben Bolker]

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