莫兰指数检验的代码-基于matlab的实现

clear; clc; A = csvread('2000.csv');%列数据 W1 = csvread('01jz.csv');%01邻接矩阵 w = normw(W1);%矩阵标准化 for k=1:18%时间跨度我这里用的是2001年至2018年 x=A(1+(k-1)*21:k*21);%1表示从第一个开始比如k=1，选用1-21单元格的数据，表示有21个地区个数 n=21;%地区个数 s=var(x,1); m=mean(x); y=0; a=0; for i=1:1:n for j=1:1:n if (i==j) continue end y=y+w(i,j)*(x(i)-m)*(x(j)-m); a=a+w(i,j); end end moran=y/(s*a); k21=0; k22=0; w0=0; w1=0; w2=0; for i=1:1:n k21=k21+(x(i)-m)^4; k22=k22+(x(i)-m)^2; for j=1:1:n w0=w0+w(i,j); w1=w1+(w(i,j)+w(j,i))^2; end end for i=1:1:n w12(i)=0; w21(i)=0; for j=1:1:n w12(i)=w12(i)+w(i,j); w21(i)=w21(i)+w(j,i); end w2=w2+(w12(i)+w21(i))^2; end k2=(n*k21)/(k22^2); w1=w1/2; ei=-1/(n-1); vari=(n*((n^2+3*n+3)*w1-n*w2+3*w0^2)-k2*((n^2-n)*w1-2*n*w2+6*w0^2))/(w0^2*(n-1)*(n-2)*(n-3))-ei^2; z=(moran-ei)/vari^(1/2); P = 2 * (1-normcdf(abs(z), 0, 1)); I(:,k)=[moran z P] end