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) 中的这些点生成条件模拟，然后简单地依赖以该点为条件的平均预测美国能源部。

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

Kriging / Gaussian Process Conditional Simulations in Matlab

matlab

gaussian

kriging