Gaussian process realization + measurements using interp2

Realization of a stochastic process is often called a sample path.

 (http://en.wikipedia.org/wiki/Sample-continuous_process)

Following is a wiki-definition:

Let (Ω, Σ, P) be a probability space. Let X : I × Ω → S be a stochastic process, where the index set I and state space S are both topological spaces. Then the process X is called sample-continuous (or almost surely continuous, or simply continuous) if the map X(ω) : I → S is continuous as a function of topological spaces for P-almost all ω in Ω.

In many examples, the index set I is an interval of time, [0, T] or [0, +∞), and the state space S is the real line or n-dimensionalEuclidean space Rⁿ.

Though it seems fairly complicated concept, in Gaussian processes, its sample path can easily be found or implemented using several Matlab lines. In this post, I'll focus on implementing GP realization in 2-dimensional space. Extending to higher dimension can easily be done. Following is the resulting figure of GP realization. 

Furthermore, implementing measuring sensory values from the GP-realized field using intperp2 function in Matlab will be introduced. This realization can be used as a reference sensory field in sensor network problems, data fusion problems, and so forth. 

1. main

%%

clc; clear; close all;

 

% 1. Make field using GP realization

stage = init_stage([0 10 0 10], [30 30]);

hyp.sig2f = 1;

hyp.sig2x = 1;

hyp.sig2w = 1E-2;

kernel = @(x, y)(kernel_se(x, y, hyp));

field = gp_realization(stage, kernel);

 

% 2. Measure from field 

nr_estimate = 50;

x_rand = stage.xmin + (stage.xmax-stage.xmin)*rand(nr_estimate, 1);

y_rand = stage.ymin + (stage.ymax-stage.ymin)*rand(nr_estimate, 1);

z_rand = zeros(nr_estimate, 1);

for zi = 1:nr_estimate

    z_rand(zi) = interp2(stage.xgrid, stage.ygrid, field.zdata, ...

            x_rand(zi), y_rand(zi) );

end

fprintf('data @(%.3f, %.3f) \t=> %.3f \n', x_rand, y_rand, z_rand);

 

% 3. Plot 

figure();pcolor(stage.xgrid, stage.ygrid, field.zdata);

set(gca,'YDir','normal');  % Matrix and Graph has different y-direction 

colormap jet; colorbar; caxis(field.zrange); hold on; 

shading interp; lighting phong; 

for ei = 1:nr_estimate

    plot(x_rand(ei), y_rand(ei), 'o', 'LineWidth',2,...

        'MarkerEdgeColor','k',...

        'MarkerSize',10);

    str = sprintf('   %.3f', z_rand(ei)); 

    text(x_rand(ei), y_rand(ei), str);

end

hold off; title('Gaussain Process Realization and Measurements', 'FontSize', 15);

2. init_stage

function stage = init_stage(axis_info, res_info)

%

% [input]

% axis_info = [xmin xmax ymin ymax];

% res_info = [xres, yres]

%

%

% [output]

% stage.xmin

% stage.xmax

% stage.ymin

% stage.ymax

% stage.xres

% stage.yres

% stage.xgrid

% stage.ygrid

% stage.points

% stage.nr_points

%

stage.xmin = axis_info(1);

stage.xmax = axis_info(2);

stage.ymin = axis_info(3);

stage.ymax = axis_info(4);

 

stage.xres = res_info(1);

stage.yres = res_info(2);

stage.xgrid = linspace(stage.xmin, stage.xmax, stage.xres);

stage.ygrid = linspace(stage.ymin, stage.ymax, stage.yres);

[mesh_X, mesh_Y] = meshgrid(stage.xgrid, stage.ygrid); 

 

% (x, y) points in column direction 

stage.points = [mesh_X(:) mesh_Y(:)]; 

stage.nr_points = length(stage.points);

3. gp_realization

function field = gp_realization(stage, kernel)

%

% [input]

% stage

% kernel

%

% [output]

% field

% 

K = kernel(stage.points, stage.points);

z = chol(K)'*randn(stage.nr_points, 1);

Z = reshape(z, stage.yres, stage.xres);

 

field.xmin = stage.xmin;

field.xmax = stage.xmax;

field.ymin = stage.ymin;

field.ymax = stage.ymax;

 

field.xgrid = stage.xgrid;

field.ygrid = stage.ygrid;

field.zdata = Z - mean(mean(Z));

field.zrange = [min(min(Z)) max(max(Z))];

4. kernel_se

function K = kernel_se(X1, X2, hyp)

 

nr_x1 = size(X1, 1);

nr_x2 = size(X2, 1);

dim = size(X1, 2);

 

% Sending error is two input data have different dimensions

assert(dim == size(X2, 2), 'Dimension of X1 and X2 should be same');

 

% Initialize kernel matrix (or vector)

K = zeros(nr_x1, nr_x2);

 

for idx1 = 1:nr_x1

    for idx2 = 1:nr_x2

        x1 = X1(idx1, :);

        x2 = X2(idx2, :);

        dist = norm(x1-x2);

        kval = hyp.sig2f*exp(-(dist^2)/hyp.sig2x);

        if (nr_x1==nr_x2) && isequal(x1, x2)

            kval = kval + hyp.sig2w;

        end

        K(idx1, idx2) = kval;

    end

end