Mikael Anderson

2006-May-11 16:04 UTC

### [R] 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]]

Ben Bolker

2006-May-13 12:07 UTC

### [R] Simulating scalar-valued stationary Gaussian processes

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