python求解全局莫兰指数和局部莫兰指数

时间：2023-09-28 05:07:02 18s弯头连接器

python解决全局莫兰指数和局部莫兰指数

1 数据简介

类别	反距离矩阵文件	属性值文件
名称	adj.csv	attribute.csv
规模	520*520	1*520
说明	无标题行和列	无标题行和列

2 相关代码

# -*- coding: utf-8 -*- import matplotlib.pyplot as plt import numpy as np import os import pandas as pd   def moranI(W, X):     ''' W:空间权重矩阵 X:观测值矩阵 归一化空间权重矩阵后moran检验，实例https://bbs.pinggu.org/thread-3568074-1-1.html '''     #W=W.astype(float)     W = np.array(W)     X = np.array(X)     X = X.reshape(1, -1)     W = W / W.sum(axis=1)  # 归一化     n = W.shape[0]  # 空间单元数     Z = X - X.mean()  # 离差阵     S0 = W.sum()     S1 = 0     for i in range(n):         for j in range(n):             S1  = 0.5 * (W[i, j]   W[j, i]) ** 2     S2 = 0     for i in range(n):         2 += (W[i, :].sum() + W[:, i].sum()) ** 2
    # 全局moran指数
    I = np.dot(Z, W)
    I = np.dot(I, Z.T)
    I = n / S0 * I / np.dot(Z, Z.T)
    # 在正太分布假设下的检验数
    EI_N = -1 / (n - 1)
    VARI_N = (n ** 2 * S1 - n * S2 + 3 * S0 ** 2) / (S0 ** 2 * (n ** 2 - 1)) - EI_N ** 2
    ZI_N = (I - EI_N) / (VARI_N ** 0.5)
    # 在随机分布假设下检验数
    EI_R = -1 / (n - 1)
    b2 = 0
    for i in range(n):
        b2 += n * Z[0, i] ** 4
    b2 = b2 / ((Z * Z).sum() ** 2)
    VARI_R = n * ((n ** 2 - 3 * n + 3) * S1 - n * S2 + 3 * S0 ** 2) - b2 * (
                (n ** 2 - n) * S1 - 2 * n * S2 + 6 * S0 ** 2)
    VARI_R = VARI_R / (S0 ** 2 * (n - 1) * (n - 2) * (n - 3)) - EI_R ** 2
    ZI_R = (I - EI_R) / (VARI_R ** 0.5)
    # 计算局部moran指数
    Ii = list()
    for i in range(n):
        Ii_ = n * Z[0, i]
        Ii__ = 0
        for j in range(n):
            Ii__ += W[i, j] * Z[0, j]
        Ii_ = Ii_ * Ii__ / ((Z * Z).sum())
        Ii.append(Ii_)
    Ii = np.array(Ii)
    # 局部检验数
    ZIi = list()
    EIi = Ii.mean()
    VARIi = Ii.var()
    for i in range(n):
        ZIi_ = (Ii[i] - EIi) / (VARIi ** 0.5)
        ZIi.append(ZIi_)
    ZIi = np.array(ZIi)
    ''' # moran散点图 # 用来正常显示中文标签 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示负号 plt.rcParams['axes.unicode_minus'] = False fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.spines['top'].set_color('none') ax.spines['right'].set_color('none') ax.xaxis.set_ticks_position('bottom') ax.spines['bottom'].set_position(('data', 0)) ax.yaxis.set_ticks_position('left') ax.spines['left'].set_position(('data', 0)) WZ = np.dot(Z, W) ax.scatter(Z, WZ, c='k') x1 = range(int(Z.min()), int(Z.max() + 1)) y1 = range(int(Z.min()), int(Z.max() + 1)) ax.plot(x1, y1, 'k--', label='x=y') x2 = list(range(int(Z.min()), int(Z.max() + 1))) y2 = np.array(x2) * I[0][0] ax.plot(x2, y2, 'k-', label='I*x=y') ax.legend(loc='upper right') imgPath = os.path.join(os.getcwd(), '莫兰散点图.png') ax.set_title('莫兰散点图') fig.savefig(imgPath) '''
    return { 
        
        'I': { 
        'value': I[0, 0], 'desc': '全局moran指数'},
        'ZI_N': { 
        'value': ZI_N[0, 0], 'desc': '正太分布假设下检验数'},
        'ZI_R': { 
        'value': ZI_R[0, 0], 'desc': '随机分布假设下检验数'},
        'Ii': { 
        'value': Ii, 'desc': '局部moran指数'},
        'ZIi': { 
        'value': ZIi, 'desc': '局部检验数'},
        #'img': {'path': imgPath, 'desc': '莫兰散点图路径'}
    }


if __name__ == "__main__":

    df = pd.read_csv('adj.csv',header=None,index_col=None)
    w=np.array(df)
    df1=pd.read_csv('attribute.csv',header=None,index_col=None)
    x = np.array(df1)
    print(len(data),len(x[0]))
    result=moranI(w, x)
    print(result['I'],result['ZI_N'])

运行结果如下：

{ 
        'value': 0.18758732033404646, 'desc': '全局moran指数'} { 
        'value': 31.11189591476361, 'desc': '正太分布假设下检验数'}

3 matlab相关方法

matlab求解全局莫兰指数：
①文件名为moransi.m
②直接调用传入参数即可



function global_moran_test(x0,w) 

row = size(x0,2);
moran.mean = zeros(row,1);
moran.num = zeros(row,1);
moran.stdev = zeros(row,1);
moran.index = zeros(row,1);
moran.z = zeros(row,1);
moran.p = zeros(row,1);

moran.globalresult = zeros(row,3);

for r = 1 : 1 : row
    x = x0(:,r);
    moran.mean(r) = mean(x);
    moran.num(r) = size(x,1);
    moran.stdev(r) = std(x);

    z_x = (x - moran.mean(r)) / moran.stdev(r);
    sum_wij = 0;
    s = 0;
    s1 = 0;
    s2 = 0;
    m2 = 0;
    m4 = 0;
    for a = 1 : 1 : moran.num(r)
        w_i = 0;
        w_j = 0;
        m2 = m2 + (x(a,1) - moran.mean(r))^2;
        m4 = m4 + (x(a,1) - moran.mean(r))^4;
        for b = 1 : 1 : moran.num(r)
            sum_wij = sum_wij + (w(a,b) * z_x(a,1) * z_x(b,1));
            s = s + w(a,b);
            s1 = s1 + (w(a,b) + w(b,a))^2;
            w_i = w_i + w(a,b);
            w_j = w_j + w(b,a);
        end
        s2 = s2 + (w_i + w_j)^2;
    end
    m2 = m2 / moran.num(r);
    m4 = m4 / moran.num(r);
    b2 = m4 / (m2^2);
    sum_i2 = 0;
    for a = 1 : 1 : moran.num(r)
        sum_i2 = sum_i2 + (z_x(a,1) * z_x(a,1));
    end
    moran.index(r) = (moran.num(r) * sum_wij) / (s * sum_i2);
    
    n = moran.num(r);
    temp_1 = n * (n^2 - 3 * n + 3) * s1 - (n * s2) + (3 * s^2);
    temp_2 = b2 * ((n^2 - n) * s1 - (2 * n * s2) + (6 * s^2));
    den = (n - 1) * (n - 2) * (n - 3) * (s^2);
    sd = (temp_1 - temp_2) / den - (1 / (n-1))^2;
    sd = sqrt(sd);

    e = -1 / (n - 1);

    moran.z(r) = (moran.index(r) - e) / sd;

    moran.p(r) = 1 - normcdf(moran.z(r));
end

moran.globalresult(:,1) = moran.index;
moran.globalresult(:,2) = moran.z;
moran.globalresult(:,3) = moran.p;

fprintf('%6s %12s %18s %24s\r\n','t','Moran','z','p');
for i = 1 : 1 : row
    fprintf('%6.3f %12.3f %18.3f %24.3\n',i,moran.index(i),moran.z(i),moran.p(i));
end
end

锐单商城拥有海量元器件数据手册、IC替代型号，打造电子元器件IC百科大全！

python求解全局莫兰指数和局部莫兰指数

相关文章