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

2 Answers 2

%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

Hello There is a video for mathmonk in youtube which explain how to do that: https://www.youtube.com/watch?v=e7as_wY0hUI&list=HL1361329139&feature=mh_lolz Regards

share|improve this answer
    
Please pay your attention to the part starting from "Also, ..." in my question. This is exactly what Mathematical Monk is showing in the video. I would like to do the same, but using the decomposition above. –  Ivan Feb 20 '13 at 6:47
1  
Please try to describe as much here as possible in order to make the answer self-contained. Links are fine as support, but they can go stale and then an answer which is nothing more than a link loses its value. –  robjohn Feb 20 '13 at 9:21
    
I thought that the Mathematical Monk used this composition, as I mentioned above. He obtained the eigenvalues and vectors and did the sampling for finite value of n. –  ibrahim Feb 20 '13 at 10:25
    
%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); 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); –  ibrahim Feb 22 '13 at 0:20

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.