“R”的“波浪号”模型公式应用于Matlab？答案

【问题标题】："Tilde" model formula of "R" applied to Matlab?“R”的“波浪号”模型公式应用于Matlab？
【发布时间】：2020-08-07 16:43:31
【问题描述】：

我现在面临一个困扰我很多的问题...我正在尝试使用 Matlab 和“BMElib2.0”工具箱从著名的“Meuse”数据（包“sp”中“R”）。问题是“lscov”函数返回了一些错误（我认为是这样）并且我的残差被 Matlab 高估了（见附件https://i.stack.imgur.com/8hXc8.png）。它在 y 轴上超过 0.5！但是，我应该得到这个（见附件https://i.stack.imgur.com/v0bU8.png）。这是 2007 Minasny & McBratney 文章“使用带有 Matérn 协方差函数的 EBLUP 进行土壤特性的空间预测”的标题，如您所见，残差不超过 0.3。我试图从“R”中得到结果，它与 Minasny & McBratney 相同（见附件https://i.stack.imgur.com/0a4YZ.png）。蓝色圆圈（残差数据）在 y 轴上的峰值为 0.3。我发现“R”使用了一个名为“模型公式”的特殊命令，它允许为另一个变量指定一个回归量。它是下面代码中的“波浪号”。我需要在 log(zinc) 和 sqrt(dist) 之间建立依赖关系。

我的问题很简单：Matlab 是否有一个“模型公式”-ish 命令以获得与“R”和 Minasny & McBratney 文章相同的残差？ 或者也许是我的“lscov”功能不适合这个？

提前感谢您提供的任何帮助！

我的“R”代码（我感兴趣的“波浪号”在“v1”命令中可见）：

library(lattice)
library(gstat)
library(sp)
load(system.file("data", "meuse.rda", package = "sp"))
v1 = variogram(log(zinc) ~ sqrt(dist), locations = ~x + y, data = meuse, width=50, cutoff=2000)
v2 = variogram(log(zinc) ~ 1, locations = ~x + y, data = meuse, width=50, cutoff=2000)
m1 = fit.variogram(v1, vgm(psill = 0.1089, "Mat", range = 40, nugget = 0.084, kappa = 8)) 
m2 = fit.variogram(v2, vgm(psill = 4.9335, "Exp", range = 1412, nugget = 0.086, kappa = 1)) 
plot(gamma~dist, v2, ylim = c(0, 1.05*max(v2$gamma)),col='red', ylab = 
       'semivariance', xlab = 'distance') 
lines(variogramLine(m2, 2000), col='red',ylab ='',xlab='') 
points(gamma~dist, v1, col='blue') 
lines(variogramLine(m1, 2000), col='blue')

我的 Matlab 代码用于残差（“BMElib2.0”用于获取变异函数）：

% "coord" (coordinates of the 155 points) and "z" (log of zinc concentration for these 155 points) 
X_Zn = [ones(size(coord(:,2))) coord(:,2) coord(:,1)];
b_derive_Zn = lscov(X_Zn,z);
Zn_derive = X_Zn*b_derive_Zn;
Zn_vec = z - Zn_derive;
figure
plZn_vec = (0:40:2000);
[dZn_vec,vZn_vec,oZn_vec] = vario(coord,Zn_vec,plZn_vec,'kron');
plot(dZn_vec,vZn_vec,'r*');
xlabel('Distance (m)','FontSize',11); ylabel('Variogram','FontSize',11);

【问题讨论】：

标签： tilde

【解决方案1】：

所以，我在回答自己：“lscov”非常适合，但我的 Matlab 代码是错误的。我试过了：

X_Zn = [ones(size(coord(:,2))) coord(:,2) coord(:,1)];

但这是不正确的。坐标在这里并不重要，重要的是“sqrt(ndist)”！所以，正确的代码是：

X_Zn = [ones(size(coord(:,2))) sqrt(ndist)];

运行脚本会给出正确的残差变异函数 :-)...希望它可以帮助其他用户！

【讨论】：