Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Assume $X(\omega, t) \sim \mathcal{N}(0, K(\cdot, \cdot))$ is a real-valued, centered Gaussian process on $R^2$, i.e., $X: \Omega \times R^2 \to R$. Let the covariance function of the process be stationary, i.e., $K(t, s) = K(t - s)$, $\forall t, s \in R^2$. More specifically, let it be exponential: $K(t, s) = \exp(-\|t - s\|/L)$ where $\|\cdot\|$ stands for the Euclidean distance, and $L$ denotes the correlation length. Using the Karhunen-Loève expansion, $X$ can be represented as the following summation:

$X(t, \omega) = \sum_{i=1}^\infty \sqrt{\lambda_i} \phi_i(t) \xi_i(\omega)$

where $\xi_i$ are independent, standard Gaussians, and $\phi_i$ and $\lambda_i$ are the eigenfunctions and eigenvalues of the covariance function $K$, which, in our case, can be computed analytically. For practical computations, the expansion above is truncated to keep, say, $n$ terms.

The question is: How to draw samples of the process on $R^2$, i.e., on a two-dimensional grid of fixed points, using the first $n$ eigenfunctions and eigenvalues of $K$? The problem is that the eigenfunctions are somewhat one-dimensional, and I cannot figure out how to properly use them to get realizations of $X$ on a plane. To be honest, the analytical solution that I claimed previously was initially obtained for processes on $R$, so I am not even sure that it applies directly to $R^2$. When I plug in the norms of points' coordinates into $\phi_i$, I get radial patterns, which is presumably not correct.

Also, it is easy to simulate $X$ by constructing the corresponding covariance matrix for a set of points, which boils down to drawing samples from a multivariate Gaussian distribution. I am not interested in this approach.

I would be very grateful for any ideas. Thank you.

share|improve this question

1 Answer 1

%choose a kernel (covariance function)

kernel = 2;
switch kernel
    case 1; k = @(x,y) 1*x'*y; %linear
    case 2; k = @(x,y) exp(-100*(x-y)'*(x-y));
    case 3; k = @(x,y) exp(-1*sqrt((x-y)'*(x-y)));
end

%choose points at which to sample
points = (0:0.02:1)';
[U,V] = meshgrid(points, points);
x=[U(:) V(:)]';
n=size(x,2);

%construct the covariance matrix
C = zeros(n,n);
for i = 1:n
    for j = 1:n
        C(i,j) = k(x(:,i), x(:,j));
    end
end
%sample from the Gaussian process at these points
u = randn(n,1);
[A,S,B]= svd(C);
z = A*sqrt(S)*u;

%plot
figure(2); clf;
Z = reshape(z, sqrt(n), sqrt(n));
surf(U,V,Z);
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.