Условное моделирование процесса Кригинга / Гаусса в Matlab

Question

Я хотел бы выполнить условное моделирование для моделей гауссовского процесса (GP) в Matlab. Я нашел учебное пособие Мартина Коларжа ( 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)

Учебное пособие генерирует условное моделирование с использованием метода прямого моделирования, основанного на разложении ковариационной матрицы. Насколько я понимаю, есть несколько методов генерации условных симуляций, которые могут быть лучше, когда число точек симуляции велико, таких как кондиционирование Кригингом с использованием локальной окрестности. Я нашел информацию о нескольких методах в Ж.-П. Чилес и П. Дельфинер, "Глава 7 - Условное моделирование", в Geostatistics: моделирование пространственной неопределенности, второе издание, John Wiley & Sons, Inc., 2012, стр. 478–628.

Существует ли существующий набор инструментов Matlab, который можно использовать для условного моделирования? Мне известны DACE, GPML и mGstat ( http://mgstat.sourceforge.net/). Я считаю, что только mGstat предлагает возможность выполнять условное моделирование. Тем не менее, mGstat также, кажется, ограничивается только 3D-моделями, и я заинтересован в моделях более высокого измерения.

Кто-нибудь может посоветовать, как начать выполнять условное моделирование с помощью существующего набора инструментов, такого как GPML?

================================================== ================= РЕДАКТИРОВАТЬ

Я нашел еще несколько наборов инструментов Matlab: STK, ScalaGauss, ooDACE

Похоже, STK способен к условному моделированию с использованием декомпозиции ковариационной матрицы. Тем не менее, ограничено умеренным количеством (может быть, несколько тысяч?) Точек симуляции из-за факторизации Холецкого.

Answer 1

Я использовал набор инструментов STK и рекомендую его другим:

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

Я обнаружил, что если вам нужно условное моделирование в большом количестве точек, то вы можете рассмотреть возможность создания условного моделирования в точках большого плана эксперимента (DoE), а затем просто полагаться на среднее предсказание, условное для этого DoE.