Получение теоретической вариограммы из параметризованной ковариационной функции (СТК)

Question

Я использую набор инструментов STK в течение нескольких дней для кригинга полей параметров окружающей среды, то есть в геостатистическом контексте.

Я нахожу набор инструментов очень хорошо реализованным и полезным (большое спасибо авторам!), И предсказания, которые я получаю через STK, на самом деле кажутся хорошими; однако я не могу визуализировать модель вариограммы, основанную на результатах STK (т. е. оценочные параметры для гауссовских функций процесса / ковариации).

Я прилагаю пример рисунка, показывающий эмпирическую вариограмму для простого 1D-теста и модель гауссовой вариограммы (как обычно используется в геостатистике, см. Также рисунок), подгоненную непосредственно к этим данным. На рисунке далее показана модель вариограммы, основанная на результатах STK, т.е. с использованием предварительно оцененных параметров модели (model.param от stk_param_estim) получить ковариацию K на целевой сетке расстояний запаздывания, а затем преобразовать K в полувариантность (в соответствии с известным соотношением semivar = K0-K, где K0 - ковариация при нулевом запаздывании). Я прилагаю простой сценарий для воспроизведения рисунка и детализации попытки преобразования.

Как вы можете видеть на рисунке, это не помогает. Я пробовал несколько других простых примеров и наборов данных STK, но модели, полученные с помощью STK и прямой подгонки, никогда не соглашаются, и на самом деле обычно выглядят гораздо более разными, чем в примере (т.е. диапазон часто кажется очень отличным, в дополнение к sill/sigma2; раскомментируйте строку 12 в скрипте, чтобы увидеть другой пример). Я также попытался ввести преобразованные параметры STK в геостатистическую модель (также в сценарии), однако вывод идентичен результату, основанному на преобразовании K выше.

Я был бы очень благодарен за вашу помощь!

Рисунок, иллюстрирующий несоответствие между вариограммами, основанное на прямом соответствии и преобразовании результатов STK

% Code to reproduce the figure illustrating my problem of getting
% variograms from STK output. The only external functions needed are those
% included with STK.
% TEST DATA - This is simply a monotonic part of the normal pdf
nugget = 0;
X = [0:20]'; % coordinates
% X = [0:50]'; % uncomment this line to see how strongly the models can deviate for different test cases
V =  normpdf(X./10+nugget,0,1); % observed values

covmodel = 'stk_gausscov_iso'; % covar model, part of STK toolbox
variomodel = 'stk_gausscov_iso_vario'; % variogram model, nested function

% GET STRUCTURE FOR THE SELECTED KRIGING (GAUSSIAN PROCESS) MODEL
nDim = size(X,2);
model = stk_model (covmodel, nDim);
model.lognoisevariance = NaN; % This makes STK fit nugget

% ESTIMATE THE PARAMETERS OF THE COVARIANCE FUNCTION
[param0, model.lognoisevariance] = stk_param_init (model, X, V); % Compute an initial guess for the parameters of the covariance function (param0)
model.param = stk_param_estim (model, X, V, param0); % Now model the covariance function

% EMPIRICAL SEMIVARIOGRAM (raw, binning removed for simplicity)
D = pdist(X)';
semivar_emp = 0.5.*(pdist(V)').^2;

% THEORETICAL SEMIVARIOGRAM FROM STK
% Target grid of lag distances
DT = [0:1:100]';
DT_zero = zeros(size(DT));
% Get covariance matrix on target grid using STK estimated pars
pairwise = true;
K = feval(model.covariance_type, model.param, DT, DT_zero, -1, pairwise);
% convert covariance to semivariance, i.e. G = C(0) - C(h)
sill = exp(model.param(1));
nugget = exp(model.lognoisevariance);
semivar_stk = sill - K + nugget; % --> this variable is then plotted

% TEST: FIT A GAUSSIAN VARIOGRAM MODEL DIRECTLY TO THE EMPIRICAL SEMIVARIOGRAM
f = @(par)mseval(par,D,semivar_emp,variomodel);
par0 = [10 10 0.1]; % initial guess for pars
[par,mse] = fminsearch(f, par0); % optimize
semivar_directfit = feval(variomodel, par, DT); % evaluate

% TEST 2: USE PARS FROM STK AS INPUT TO GAUSSIAN VARIOGRAM MODEL
par(1) = exp(model.param(1)); % sill, PARAM(1) = log (SIGMA ^ 2), where SIGMA is the standard deviation,
par(2) = sqrt(3)./exp(model.param(2)); % range, PARAM(2) = - log (RHO), where RHO is the range parameter. --- > RHO = exp(-PARAM(2))
par(3) = exp(model.lognoisevariance); % nugget
semivar_stkparswithvariomodel = feval(variomodel, par, DT);

% PLOT SEMIVARIOGRAM
figure(); hold on;
plot(D(:), semivar_emp(:),'.k'); % Observed variogram, raw
plot(DT, semivar_stk,'-b','LineWidth',2); % Theoretical variogram, on a grid
plot(DT, semivar_directfit,'--r','LineWidth',2); % Test direct fit variogram
plot(DT,semivar_stkparswithvariomodel,'--g','LineWidth',2); % Test direct fit variogram using pars from stk
legend('raw empirical semivariance (no binned data here for simplicity) ',...
    'Gaussian cov model from STK, i.e. exp(Sigma2) - K + exp(lognoisevar)',...
    'Gaussian semivariogram model (fitted directly to semivariance)',...
    'Gaussian semivariogram model (using transformed params from STK)');
xlabel('Lag distance','Fontweight','b');
ylabel('Semivariance','Fontweight','b');

% NESTED FUNCTIONS
% Objective function for direct fit
function [mse] = mseval(par,D,Graw,variomodel)
Gmod = feval(variomodel, par, D);
mse = mean((Gmod-Graw).^2);
end
% Gaussian semivariogram model.
function [semivar] = stk_gausscov_iso_vario(par, D) %#ok<DEFNU>
% D : lag distance, c : sill, a : range, n : nugget
c = par(1); % sill
a = par(2); % range
if length(par) > 2, n = par(3); % nugget optional
else, n = 0; end
semivar = n + c .* (1 - exp( -3.*D.^2./a.^2 )); % Model
end

Answer 1

Нет ничего плохого в том, как вы вычисляете вариограмму.

Чтобы понять фигуру, которую вы получаете, учтите, что:

1. Параметры модели оцениваются в STK с использованием (ограниченного) метода максимального правдоподобия, а не методом наименьших квадратов на полувариограмме.
2. Для очень гладких стационарных случайных полей, наблюдаемых через короткие промежутки времени, не следует ожидать, что теоретическая вариограмма будет совпадать с эмпирической вариограммой, с или без биннинга. Причина этого заключается в том, что наблюдения и, следовательно, квадратные различия в этом случае очень коррелируют.

Чтобы убедить себя во втором пункте, вы можете запустить следующий скрипт несколько раз:

% a smooth GP
model = stk_model (@stk_gausscov_iso, 1);
model.param = log ([1.0, 0.2]);  % unit variance

x_max = 20;  x_obs = x_max * rand (50, 1);

% Simulate data
z_obs = stk_generate_samplepaths (model, x_obs);

% Empirical semivariogram (raw, no binning)
h = (pdist (double (x_obs)))';
semivar_emp = 0.5 * (pdist (z_obs)') .^ 2;

% Model-based semivariogram
x1 = (0:0.01:x_max)';
x0 = zeros (size (x1));
K = feval (model.covariance_type, model.param, x0, x1, -1, true);
semivar_th = 1 - K;

% Figure
figure;  subplot (1, 2, 1);  plot (x_obs, z_obs, '.');
subplot (1, 2, 2);  plot (h(:), semivar_emp(:),'.k');  hold on;
plot (x1, semivar_th,'-b','LineWidth',2);
legend ('empirical', 'model');  xlabel ('lag');  ylabel ('semivar');

Дальнейшие вопросы по оценке параметров для моделей гауссовских процессов, вероятно, следует задавать в отношении перекрестной проверки, а не переполнения стека.