R help - May 2006 - Simulating scalar-valued stationary Gaussian processes

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.

/Mikeal

	[[alternative HTML version deleted]]

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