1

几天来，我一直在使用 STK 工具箱，用于环境参数字段的克里金法，即在地统计环境中。

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

我附上了一个示例图，显示了一个简单的 1D 测试用例的经验半变异函数和直接拟合该数据的高斯半变异函数模型（通常用于地统计学，另请参见图）。该图进一步显示了基于 STK 输出的半变异函数模型，即使用先前估计的模型参数 ( model.paramfrom 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
4

1 回答 1

0

计算半变异函数的方式没有任何问题。

要了解您获得的数字，请考虑：

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

关于高斯过程模型的参数估计的进一步问题可能应该在 Cross-Validated 而不是 Stack Overflow 上提出。

于 2018-03-03T21:10:26.717 回答