找到你要的答案

Q:Kriging / Gaussian Process Conditional Simulations in Matlab

Q:克里格/高斯过程的条件模拟在Matlab

I would like to perform conditional simulations for Gaussian process (GP) models in Matlab. I have found a tutorial by 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)

The tutorial generates the conditional simulations using a direct simulation method based on covariance matrix decomposition. It is my understanding that there are several methods of generating conditional simulations that may be better when the number of simulation points is large such as conditioning by Kriging using a local neighborhood. I have found information regarding several methods in J.-P. Chilès and P. Delfiner, “Chapter 7 - Conditional Simulations,” in Geostatistics: Modeling Spatial Uncertainty, Second Edition, John Wiley & Sons, Inc., 2012, pp. 478–628.

Is there an existing Matlab toolbox that can be used for conditional simulations? I am aware of DACE, GPML, and mGstat (http://mgstat.sourceforge.net/). I believe only mGstat offers the capability to perform conditional simulations. However, mGstat also seems to be limited to only 3D models and I am interested in higher dimensional models.

Can anybody offer any advice on getting started performing conditional simulations with an existing toolbox such as GPML?

=================================================================== EDIT

I have found a few more Matlab toolboxes: STK, ScalaGauss, ooDACE

It appears STK is capable of conditional simulations using covariance matrix decomposition. However, is limited to a moderate number (maybe a few thousand?) of simulation points due to the Cholesky factorization.

我想执行条件模拟高斯过程(GP)模型在Matlab。我已经通过马丁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.孩子è和P. Delfiner的几种方法发现的信息,“7章条件模拟,“Geostatistics:空间的不确定性,第二版造型,John Wiley &;子公司,2012,页478–628。

有一个现有的MATLAB工具箱,可用于条件模拟?我知道鲮鱼,吸收边界条件,并mgstat(HTTP:/ / mgstat。SourceForge。网/)。我相信只有mgstat提供执行条件的模拟能力。然而,mgstat似乎也仅限于3D模型和我在高维模型感兴趣。

谁能提供任何意见,开始与现有的工具箱如GPML执行条件模拟?

=================================================================== EDIT

我发现几个MATLAB工具箱:STK、ScalaGauss、oodace

看来STK可以使用协方差矩阵分解的条件模拟。然而,仅限于适度的数字(也许几千?)模拟点由于Cholesky分解。

answer1: 回答1:

I used the STK toolbox and I recommend it for others:

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

I found that if you need conditional simulations at a large number of points then you might consider generating a conditional simulation at the points in a large design of experiment (DoE) and then simply relying on the mean prediction conditional on that DoE.

我用STK工具箱和我推荐它的人:

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

我发现,如果你需要有条件的模拟在大量的点,那么你可能会考虑产生一个条件模拟的点在一个大的设计实验(DOE),然后简单地依赖于平均预测条件的DOE。

matlab  gaussian  kriging