用matlab求moran并作图的具体步骤

现在有从论坛上下载的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

这个代码确实只能做单变量的哈，我也是这几天在研究，代码有一些小问题，个人觉得

请问你选择解决了吗，同问题

同问，你的问题解决了没