2

我想在 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 分解，模拟点的数量仅限于中等数量（可能是几千？）。

4

1 回答 1

2

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

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

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

于 2016-09-23T16:14:27.373 回答