现在有从论坛上下载的moran 的代码,如下:求问这个代码如何修正,运行总是有问题。这个代码只能做单变量的moran 吗?

function[I,I_standard,Z_I,Z_alpha,result]=moran(X,W,alpha)
%%X为列向量,W为权重矩阵
%alpha为显著性水平
n = length(X);
Xmean = mean(X);
XminusX_mean = X - Xmean;
S0 = sum(sum(W));
%一、计算,moran'I值
%1、W未标准化,为0-1矩阵
I =XminusX_mean'*W*XminusX_mean*n/(XminusX_mean'*XminusX_mean*S0);
%2、W标准化,为(0,1)矩阵
W_rowSum = sum(W');
W_rowSum_temp = W_rowSum'*ones(1,n);
W_standard = W./W_rowSum_temp;
I_standard =XminusX_mean'*W_standard*XminusX_mean/(XminusX_mean'*XminusX_mean);
%二、moran'I值的显著性检验
Z_I = moran_test(I_standard,W_standard);
Z_alpha = norminv(1-alpha,0,1);
%三、画moran散点图
%y = ax + b;线性拟合
%1、X标准化,W标准化
Xstd = std(X,1);%X的标准差
X_std = (X-Xmean)/Xstd;%X标准化
W_standard_X_standard = W_standard*X_std;
%result = myls(X_std,WX_standard);
result = moranScatterPlot(X_std, W_standard_X_standard );
title('X标准化,W标准化');
% %2、X标准化,W未标准化
% Xstd = std(X,1);%X的标准差
% X_std = (X-Xmean)/Xstd;%X标准化,X_std代表横坐标
% W_X_standard = W*X_std;%WX_standard代表纵坐标
% % result = myls(X_std,WX_standard);
% result = moranScatterPlot(X_std,W_X_standard);
% title('X标准化,W未标准化');
%3、X为未标准化,W标准化
% W_standard_X = W_standard*X;%WX_standard代表纵坐标
% % result = myls(X_std,WX_standard);
% result = moranScatterPlot(X,W_standard_X)
% title('X为未标准化,W标准化');
% % 4、X未标准化,W未标准化
% WX = W*X;%WX_standard代表纵坐标
% result = moranScatterPlot(X,WX);
% title('X未标准化,W未标准化');
Moran I 检验
function [Z_I,Z_alpha] =moran_test(I,W,alpha)
%%I为计算出来的moran'I值
%W为计算相应moran'I值的权重矩阵
n = size(W,1);
E_I = -1/(n-1);
S0 = sum(sum(W));
S1 = sum( sum( (W+W').^2 ) )/2;
S2 =sum( sum((W+W')').^2 );
Var_I =(n*n*S1-n*S2+3*S0*S0)/((n*n-1)*S0*S0)-E_I*E_I;
Z_I = (I-E_I)/sqrt(Var_I);
Z_alpha = norminv(1-alpha,0,1);
Moran I plot
function result = moranScatterPlot(X,WX)
%%X为列向量,W为权重矩阵
result = myls(X,WX);
X_lift = min(X) - 0.5*abs(min(X));
X_right = max(X) + 0.5*max(X);
xx = linspace(X_lift,X_right,1000);
yy = result(1)*xx + result(2);
plot(X,WX,'ro');
hold on
plot(xx,yy)
grid on
axis equal
end