matlab - Matlab 中的克里金/高斯过程条件模拟

Question

我想在 Matlab 中对高斯过程 (GP) 模型进行条件模拟。我找到了 Martin Kolář 的教程 ( 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)

本教程使用基于协方差矩阵分解的直接模拟方法生成条件模拟。据我了解，有几种生成条件模拟的方法，当模拟点的数量很大时，这些方法可能会更好，例如使用局部邻域进行克里金法调节。我在 J.-P. 中找到了有关几种方法的信息。 Chilès 和 P. Delfiner，“第 7 章 - 条件模拟”，《地统计学:空间不确定性建模》，第二版，John Wiley & Sons, Inc.，2012 年，第 478–628 页。

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

任何人都可以提供有关如何开始使用 GPML 等现有工具箱执行条件模拟的建议吗？

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

我又找到了几个Matlab工具箱:STK、ScalaGauss、ooDACE

看来 STK 能够使用协方差矩阵分解进行条件模拟。然而，由于 Cholesky 分解，模拟点的数量有限(也许几千？)。

Answer 1

我使用了 STK 工具箱，并向其他人推荐它:

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

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