基于Matlab的克里格/高斯过程条件模拟

Question

我想在Matlab中对高斯过程(GP)模型进行条件模拟。我找到了马丁·科拉( Martinář，http://mrmartin.net/?p=223)的教程。

sigma_f = 1.1251; %parameter of the squared exponential kernel
l =  0.90441; %parameter of the squared exponential kernel
kernel_function = @(x,x2) sigma_f^2*exp((x-x2)^2/(-2*l^2));

%This is one of many popular kernel functions, the squared exponential
%kernel. It favors smooth functions. (Here, it is defined here as an anonymous
%function handle)

% we can also define an error function, which models the observation noise
sigma_n = 0.1; %known noise on observed data
error_function = @(x,x2) sigma_n^2*(x==x2); 
%this is just iid gaussian noise with mean 0 and variance sigma_n^2s

%kernel functions can be added together. Here, we add the error kernel to
%the squared exponential kernel)
k = @(x,x2) kernel_function(x,x2)+error_function(x,x2);

X_o = [-1.5 -1 -0.75 -0.4 -0.3 0]';
Y_o = [-1.6 -1.3 -0.5 0 0.3 0.6]';

prediction_x=-2:0.01:1;

K = zeros(length(X_o));
for i=1:length(X_o)
    for j=1:length(X_o)
        K(i,j)=k(X_o(i),X_o(j));
    end
end

%% Demo #5.2 Sample from the Gaussian Process posterior
clearvars -except k prediction_x K X_o Y_o

%We can also sample from this posterior, the same way as we sampled before:
K_ss=zeros(length(prediction_x),length(prediction_x));
for i=1:length(prediction_x)
    for j=i:length(prediction_x)%We only calculate the top half of the matrix. This an unnecessary speedup trick
        K_ss(i,j)=k(prediction_x(i),prediction_x(j));
    end
end
K_ss=K_ss+triu(K_ss,1)'; % We can use the upper half of the matrix and copy it to the

K_s=zeros(length(prediction_x),length(X_o));
for i=1:length(prediction_x)
    for j=1:length(X_o)
        K_s(i,j)=k(prediction_x(i),X_o(j));
    end
end

[V,D]=eig(K_ss-K_s/K*K_s');
A=real(V*(D.^(1/2)));

for i=1:7
    standard_random_vector = randn(length(A),1);
    gaussian_process_sample(:,i) = A * standard_random_vector+K_s/K*Y_o;
end
hold on
plot(prediction_x,real(gaussian_process_sample))

set(plot(X_o,Y_o,'r.'),'MarkerSize',20)

本教程使用基于协方差矩阵分解的直接模拟方法生成条件模拟。据我所知，有几种产生条件模拟的方法，当模拟点的数目很大时，可能会更好，例如使用局部邻域的Kriging调节。我在J.-P. Chilès和P. Delfiner的“第7章-条件模拟”中发现了关于几种方法的信息，载于Geostatistics:模拟空间不确定性模型，第二版，John Wiley & Sons，Inc.，2012年，第478至628页。

现有的Matlab工具箱可以用于条件模拟吗？I知道DACE、GPML和mGstat (http://mgstat.sourceforge.net/)。我相信只有mGstat提供了执行条件模拟的能力。然而，mGstat似乎也仅限于三维模型，我对高维模型感兴趣。

是否有人对开始使用现有工具箱(如GPML? )执行条件模拟提供任何建议？

===================================================================编辑

我还找到了一些Matlab工具箱: STK，ScalaGauss，ooDACE

STK可以用协方差矩阵分解进行条件模拟。然而，仅限于一个适度的数字(可能几千？)由于Cholesky因式分解而产生的模拟点。

Answer 1

我使用了STK工具箱，并推荐给其他人：

http://kriging.sourceforge.net/htmldoc/

我发现，如果需要在大量的点上进行条件模拟，那么可以考虑在大型实验设计( DoE )中的点处生成一个条件模拟，然后简单地依赖于该DoE上的平均预测。