从参数化协方差函数（STK）中获得理论半变异函数

Question

我一直在使用STK工具箱几天，用于克服环境参数字段，即在地统计背景下。

我发现工具箱已经很好地实现并且非常有用（非常感谢作者！），而我通过STK获得的克里金预测实际上看起来很好;然而，我发现自己无法基于STK输出（即高斯过程/协方差函数的估计参数）可视化半变异函数模型。

我附上一个示例图，显示了一个简单的1D测试用例的经验半变异函数和一个直接拟合到该数据的高斯半变异函数模型（通常用于地质统计学，见图）。该图进一步显示了基于STK输出的半变异函数模型，即使用先前估计的模型参数（来自model.param的stk_param_estim）来获得滞后距离的目标网格上的协方差K，然后将K转换为半方差（根据众所周知的关系半公式= K0-K其中K0是零滞后的协方差。我附上一个简单的脚本来重现图形并详细说明尝试的转换。

正如您在图中所看到的，这并不能解决问题。我已经尝试了其他一些简单的例子和​​STK数据集，但是通过STK与直接拟合获得的模型从未达成一致，实际上通常看起来与示例中的不同（即范围通常看起来非常不同，除了sill / sigma2 ;取消注释脚本中的第12行以查看另一个示例）。我还尝试将转换后的STK参数输入到地统计模型中（也在脚本中），但是，输出与基于转换K的结果相同。

我非常感谢你的帮助！

Figure illustrating the lack of agreement between semivariograms based on direct fit vs conversion of STK output

% 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');

关于高斯过程模型的参数估计的其他问题可能应该在交叉验证而不是堆栈溢出上提出。