はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
　本を読んだ上で補助的に読んでください。　

　この記事では、グローバル・モランのIについて、図を使って解説します。

【前の内容】

【他の内容】

【今回の内容】

2.3 空間的自己相関(Global Moran's I)の可視化

　空間的自己相関(spatial autocorrelation)の指標であるグローバル・モランのI(GMI・global Moran's I)を図で確認します。
　GMIの計算については「【Python】2.3：空間的自己相関(Global Moran's I)の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」、空間重み行列(spatial weight matrix)については「【Python】2.2：空間重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import geopandas as gpd
from pysal.lib import weights
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib # 日本語の表示用
from mpl_toolkits.axes_grid1 import make_axes_locatable # カラーバーのサイズ調整用
from matplotlib.animation import FuncAnimation

グローバル・モランのIの数式表現

　まずは、GMIの定義や計算を数式で確認します。

　 $N$ 個の区域に関して、 $i$ 番目の区域のデータ(変数)を $x_i$ として、 $N$ 個のデータをまとめて $N$ 次元ベクトル $\mathbf{x}$ 、空間重み行列を $(N \times N)$ の行列 $\mathbf{W}$ とします。

$\displaystyle \mathbf{x} = \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{pmatrix} ,\ \mathbf{W} = \begin{pmatrix} w_{11} & w_{12} & \cdots & w_{1N} \\ w_{21} & w_{22} & \cdots & w_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ w_{N1} & w_{N2} & \cdots & w_{NN} \end{pmatrix}$

　GMIは、次の式で定義されます。

$\displaystyle \mathrm{GMI} = \frac{1}{\sum_{i=1}^N \sum_{j=1}^N w_{ij}} \frac{ \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x}) }{ \frac{1}{N} \sum_{i=1}^N (x_i - \bar{x})^2 }$

　ベクトルや行列を用いた計算については「2.3-4：空間的自己相関(Moran's I)の計算式の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

グローバル・モランのIの実装

　続いて、GMIの計算を行う関数を作成します。

　GMIの計算を自作関数として定義します。

# GMIの計算を実装
def culc_GMI(x, W):
    
    # データ数を取得
    N = len(x)
    
    # 偏差を計算
    x_dev = x - np.mean(x)

    # GMIを計算
    GMI = N / np.sum(W) * (x_dev.T @ W @ x_dev) / (x_dev.T @ x_dev)
    return GMI.item()

　ここまでは、他の記事で確認しました。

データの読込

　可視化の前に、例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。

　「国土数値情報ダウンロードサイト」の行政区域データと地価公示データを読み込みます。

処理コード(クリックで展開)

# ファイルパスを指定
DISTRICT_PATH = 'N03-20230101_27_GML/N03-23_27_230101.shp'
LAND_PATH     = 'L01-23_27_GML/L01-23_27.shp'

# データを読込
gdf_district = gpd.read_file(DISTRICT_PATH, encoding='shift-jis') # 行政区域
gdf_land     = gpd.read_file(LAND_PATH) # 地価公示

# 対象データを取得
gdf_district = gdf_district[['N03_003', 'N03_004', 'N03_007', 'geometry']]
gdf_district.columns = ['city1', 'city2', 'd_code', 'geometry']
gdf_land = gdf_land[['L01_024', 'L01_022', 'L01_006', 'geometry']]
gdf_land.columns = ['address', 'd_code', 'price', 'geometry']

# データを整形
dst_proj = 6668
gdf_district = gdf_district.to_crs(epsg=dst_proj)
gdf_land     = gdf_land.to_crs(epsg=dst_proj) # 空間座標系を再設定
gdf_district['city2_label'] = gdf_district['city2'].str.replace('大阪市|堺市', '', regex=True) # 地名が重なる対策用

# データを加工
gdf_land['log_price'] = np.log(gdf_land['price'])

　公示価格の対数をとった列を作成します。その他にも扱いやすいように整形しておきます。

　行政区域データと地価公示データをまとめます。

# データを統合
gdf_target = gdf_district.merge(
    gdf_land.groupby('d_code')['log_price'].mean(), 
    on='d_code'
)
gdf_target = gdf_target.dissolve(by=['d_code'], as_index=False)
gdf_target['centers'] = gdf_target['geometry'].centroid

# データフレームを整形
gdf_target = gdf_target.reindex(columns=['city1', 'city2', 'd_code', 'geometry', 'centers', 'log_price', 'city2_label'])
gdf_target

	city1	city2	d_code	geometry	centers	log_price	city2_label
0	大阪市	大阪市都島区	27102	POLYGON ((135.51481 34.72057, 135.51489 34.720...	POINT (135.52721 34.71193)	12.769969	都島区
1	大阪市	大阪市福島区	27103	POLYGON ((135.46247 34.6874, 135.46253 34.6879...	POINT (135.47534 34.6939)	13.214822	福島区
2	大阪市	大阪市此花区	27104	MULTIPOLYGON (((135.36506 34.6325, 135.35885 3...	POINT (135.42012 34.668)	12.191524	此花区
3	大阪市	大阪市西区	27106	POLYGON ((135.46667 34.67284, 135.4642 34.6747...	POINT (135.48375 34.67776)	13.746399	西区
4	大阪市	大阪市港区	27107	MULTIPOLYGON (((135.43009 34.64977, 135.4301 3...	POINT (135.45084 34.66032)	12.410058	港区
...	...	...	...	...	...	...	...
67	泉南郡	田尻町	27362	MULTIPOLYGON (((135.29259 34.40094, 135.29293 ...	POINT (135.25775 34.417)	10.923984	田尻町
68	泉南郡	岬町	27366	MULTIPOLYGON (((135.18163 34.34026, 135.18187 ...	POINT (135.15203 34.30321)	9.918761	岬町
69	南河内郡	太子町	27381	POLYGON ((135.66132 34.53314, 135.66139 34.532...	POINT (135.65269 34.51765)	10.697891	太子町
70	南河内郡	河南町	27382	MULTIPOLYGON (((135.64123 34.45932, 135.64105 ...	POINT (135.64897 34.48455)	10.327751	河南町
71	南河内郡	千早赤阪村	27383	POLYGON ((135.62578 34.47662, 135.62579 34.476...	POINT (135.64812 34.43844)	10.122623	千早赤阪村

72 rows × 7 columns

　市区町村(標準地行政区域コード)ごとに、対数公示価格の平均をとった列を作成して、ポリゴンデータと結合します。また、ポリゴンデータをまとめて、中心座標を求めます。

　以上で、作図に利用する空間データを用意できました。

グローバル・モランのIの計算

　次は、空間データからGMIを計算します。

　変数を設定します。

# 変数を抽出
x = gdf_target['log_price'].values
print(x[:5].round(2))
print(x.shape)

# 区域数を取得
N = len(x)
print(N)

[12.77 13.21 12.19 13.75 12.41]
(72,)
72

　この例では、区域ごとの平均対数公示価格を変数として用います。

　空間隣接行列を作成して、空間重み行列に変換します。

# 隣接関係を作成
wr = weights.Rook.from_dataframe(gdf_target)

# 空間隣接行列を作成
adj_mat = wr.full()[0].astype(int)
print(adj_mat)
print(adj_mat.shape)

# 空間重み行列を作成
weight_mat = adj_mat / adj_mat.sum(axis=1, keepdims=True)
np.nan_to_num(weight_mat, nan=0.0) # 列要素が全て0の場合用
print(weight_mat.round(2))
print(weight_mat.shape)

[[0 0 0 ... 0 0 0]
 [0 0 1 ... 0 0 0]
 [0 1 0 ... 0 0 0]
 ...
 [0 0 0 ... 0 1 0]
 [0 0 0 ... 1 0 1]
 [0 0 0 ... 0 1 0]]
(72, 72)
[[0.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.   0.2  ... 0.   0.   0.  ]
 [0.   0.25 0.   ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.   0.33 0.  ]
 [0.   0.   0.   ... 0.33 0.   0.33]
 [0.   0.   0.   ... 0.   0.33 0.  ]]
(72, 72)

　weights ライブラリを利用して、空間データから空間隣接行列を作成して、空間重み行列に変換します。

　GMIを計算します。

# 重みを設定
W = weight_mat.copy()

# GMIを計算
GMI = culc_GMI(x, W)
print(GMI)

0.744821273246302

　自作関数を使ってGMIを計算します。

グローバル・モランのIの図

　続いて、GMIの定義や性質を可視化します。コロプレス図については本を参照してください。
　変数と偏差の関係については「【Python】2.3：空間ラグの可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　偏差のコロプレス図とヒートマップを作成します。

# 偏差を計算
x_dev = x - np.mean(x)

# 偏差を格納
gdf_target['deviation'] = x_dev

# 色の調整用
x_min = np.floor(x_dev.min())
x_max = np.ceil(x_dev.max())

# 偏差を作図
fig, axes = plt.subplots(nrows=1, ncols=2, constrained_layout=True, 
                       figsize=(12, 10), width_ratios=[8, 2], dpi=100, facecolor='white')
fig.suptitle('log price (deviation)', fontsize=20)

# コロプレスマップを作図
ax = axes[0]
gdf_district.boundary.plot(ax=ax, linewidth=0.5, edgecolor='white') # 行政区界
gdf_target.plot(ax=ax, column='deviation', cmap='jet', vmin=x_min, vmax=x_max) # 各区域の変数
for label_x, label_y, label_str in zip(gdf_target.centers.x, gdf_target.centers.y, gdf_target.city2_label):
    ax.annotate(text=label_str, xy=(label_x-0.015, label_y-0.005), size=9) # 区域名
ax.set_xlabel('longitude')
ax.set_ylabel('latitude')
ax.set_title(f'Osaka: $N = {N}, \\bar{{x}} = {np.mean(x):.1f}, GMI = {GMI:.3f}$', loc='left')
ax.grid()
ax.set_aspect('equal', adjustable='box')

# ヒートマップを作図
ax = axes[1]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='200%', pad=0.25)
c = ax.pcolor(x_dev[:, np.newaxis], cmap='jet', vmin=x_min, vmax=x_max) # ベクトル
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=gdf_target.city2)
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('city ( $i$ )')
ax.grid()
fig.colorbar(c, cax=cax, label='$x_i - \\bar{x}$')
ax.set_aspect('equal', adjustable='box')

plt.show()

　左のコロプレス図は、偏差 $x_i - \bar{x}$ の値に応じて区域ごとを塗りつぶして示します。右のヒートマップは、全ての区域を1列に並べた図です。

　変数や隣接関係によってGMIが決まります。

　GMIの計算のヒートマップを作成します。

作図コード(クリックで展開)

# 偏差を計算
x_dev = (x - np.mean(x))[:, np.newaxis]

# 偏差の積を計算
X_dev_prod = x_dev @ x_dev.T

# 基準化済み重み付け偏差の積を計算
Z = N/np.sum(W) * W*X_dev_prod / np.sum(x_dev**2)

# 色の設定用
x_max  = np.ceil(max(x_dev.max(), abs(x_dev.min()), X_dev_prod.max(), abs(X_dev_prod.min())))
w_max  = 1.0
wx_max = max((W*X_dev_prod).max(), abs((W*X_dev_prod).min()))
z_max  = max(Z.max(), abs(Z.min()))

# GMIの計算を作図
fig, axes = plt.subplots(nrows=3, ncols=3, constrained_layout=True, 
                       figsize=(23, 21), width_ratios=[1, 11, 11], height_ratios=[10, 1, 10], 
                       dpi=100, facecolor='white')
fig.suptitle("Global Moran's I", fontsize=20)

# 偏差を作図
ax = axes[0, 0]
ax.pcolor(x_dev, cmap='jet', vmin=-x_max, vmax=x_max) # 偏差
ax.set_xticks(ticks=[0.5])
ax.set_xticklabels(labels='')
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=gdf_target.city2, size=9) # 区域名
ax.invert_yaxis() # 軸の反転
ax.set_ylabel('city ( $i$ )')
ax.set_title('$x - \\bar{x}$', fontsize=20)
ax.grid()
ax.set_aspect('equal', adjustable='box')

# 偏差を作図
ax = axes[1, 1]
ax.pcolor(x_dev.T, cmap='jet', vmin=-x_max, vmax=x_max) # 偏差
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=gdf_target.city2, size=9, rotation=90) #区域名
ax.set_yticks(ticks=[0.5])
ax.set_yticklabels(labels='')
ax.set_xlabel('city ( $j$ )')
ax.set_title('$(x - \\bar{x})^T$', fontsize=20)
ax.grid()
ax.set_aspect('equal', adjustable='box')

# 偏差の積を作図
ax = axes[0, 1]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(X_dev_prod, cmap='jet', vmin=-x_max, vmax=x_max) # 偏差の積
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title('$(x-\\bar{x}) (x-\\bar{x})^T$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$(x_i-\\bar{x}) (x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# 空間重み行列を作図
ax = axes[0, 2]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(W, cmap='seismic', vmin=-w_max, vmax=w_max) # 重み
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title('$W$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$w_{ij}$')
ax.set_aspect('equal', adjustable='box')

# 重み付け偏差の積を作図
ax = axes[2, 1]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(W*X_dev_prod, cmap='RdYlBu_r', vmin=-wx_max, vmax=wx_max) # 重み付け偏差の積
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title('$W \\odot (x-\\bar{x}) (x-\\bar{x})^T$', fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label='$w_{ij} (x_i-\\bar{x}) (x_j-\\bar{x})$')
ax.set_aspect('equal', adjustable='box')

# ラベル用の文字列を作成
fnc_label   = '$\\frac{N}{\\sum_{i=1}^N \\sum_{j=1}^N w_{ij}} '
fnc_label  += 'W \\odot (x-\\bar{x}) (x-\\bar{x})^T'
fnc_label  += '\\frac{1}{\\sum_{i=1}^N (x_i-\\bar{x})}$'
elem_label  = '$\\frac{N}{\\sum_{i=1}^N \\sum_{j=1}^N w_{ij}} '
elem_label += '\\frac{w_{ij} (x_i-\\bar{x}) (x_j-\\bar{x})}{\\sum_{i=1}^N (x_i-\\bar{x})}$'

# 基準化済み重み付け偏差の積を作図
ax = axes[2, 2]
divider = make_axes_locatable(ax)
cax = divider.append_axes('right', size='5%', pad=0.25)
c = ax.pcolor(Z, cmap='RdYlBu_r', vmin=-z_max, vmax=z_max) # 基準化済み重み付け偏差の積
ax.set_xticks(ticks=np.arange(N)+0.5)
ax.set_xticklabels(labels=np.arange(N)+1, size=6.5)
ax.set_yticks(ticks=np.arange(N)+0.5)
ax.set_yticklabels(labels=np.arange(N)+1, size=8)
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$j$')
ax.set_ylabel('$i$')
ax.set_title(fnc_label, fontsize=20)
ax.grid()
fig.colorbar(c, cax=cax, label=elem_label)
ax.set_aspect('equal', adjustable='box')

# GMIを表示
ax = axes[1, 2]
ax.text(x=0, y=0, s=f'GMI = {culc_GMI(x, W):.3}', fontsize=20)
ax.set_xlim(xmin=-1, xmax=1)

# 不要な図を非表示
axes[1, 0].axis('off')
axes[1, 2].axis('off')
axes[2, 0].axis('off')

plt.show()

　それぞれ値の意味合いが異なるのでカラーマップ(配色)を変えて、値の大小をグラデーションで示します。
　偏差ベクトル(ベクトルとスカラの差)は、正確には $\mathbf{x} - \bar{x} \mathbf{1}$ ですが、この図では(これ以上ゴチャゴチャしてもなんなので) $\mathbf{x} - \bar{x}$ と表記しています。

　左上図は、縦ベクトルと横ベクトル(転置した)の偏差と、2区域の全ての組み合わせの「偏差の積」です。偏差の積の範囲で偏差も配色しています。
　右上図は、隣接重み行列です。
　左下図は、隣接重み行列によって「重み付けした偏差の積」(左上図と右上図の要素ごとの積)です。隣接しない組み合わせは0になりGMIに影響せず、隣接数が少ない区域との組み合わせほど値が大きく(割引率が小さく)なります。
　右下図は、「隣接重み行列の総和」と「分散(偏差の2乗和の平均)」によって基準化した「重み付け偏差の積」です。総和の絶対値が1以下になります。
　右下図の総和がGMIです(空いているスペースに表示しています)。

隣接関係とGMIの関係

　簡単な例として、格子状に隣接する区域における変数とGMIの関係を図で確認します。

作図コード(クリックで展開)

　区域数を設定します。

# 一辺の区域数を指定
n = 5

# 全体の区域数を計算
N = n**2
print(N)

　格子状に並ぶ区域の一辺の区域数 n を指定して、全体の区域数 N を計算します。

　空間隣接行列を作成します。

# 空間隣接行列を作成
tilde_W = np.zeros((N, N), dtype='int')
adj_i = np.arange(n, N)
adj_j = np.arange(0, N-n)
tilde_W[adj_i, adj_j] = 1 # 上:(i, i-n)
tilde_W[adj_j, adj_i] = 1 # 下:(i, i+n)
adj_i = np.hstack([np.arange(n*i+1, n*(i+1)) for i in range(n)])
adj_j = np.hstack([np.arange(n*i, n*(i+1)-1) for i in range(n)])
tilde_W[adj_i, adj_j] = 1 # 左:(i, i-1)
tilde_W[adj_j, adj_i] = 1 # 右:(i, i+1)
print(tilde_W[:10, :10])

[[0 1 0 0 0 1 0 0 0 0]
 [1 0 1 0 0 0 1 0 0 0]
 [0 1 0 1 0 0 0 1 0 0]
 [0 0 1 0 1 0 0 0 1 0]
 [0 0 0 1 0 0 0 0 0 1]
 [1 0 0 0 0 0 1 0 0 0]
 [0 1 0 0 0 1 0 1 0 0]
 [0 0 1 0 0 0 1 0 1 0]
 [0 0 0 1 0 0 0 1 0 1]
 [0 0 0 0 1 0 0 0 1 0]]

　格子状に並ぶ区域に対応する空間隣接行列を作成します。
　一辺の区域数を $n$ とすると、区域 $i$ の上の区域番号は $i-n$ 、下の区域番号は $i+n$ 、左の区域番号は $i-1$ 、右の区域番号は $i+1$ です。ただし、区域 $i$ の場所によって、上下左右のどの区域が存在するかは異なります。

　空間重み行列を作成します。

# 隣接重み行列に変換
W = tilde_W / tilde_W.sum(axis=1, keepdims=True)
np.nan_to_num(W, 0.0) # 列要素が全て0の場合用
print(W[:10, :10].round(2))

[[0.   0.5  0.   0.   0.   0.5  0.   0.   0.   0.  ]
 [0.33 0.   0.33 0.   0.   0.   0.33 0.   0.   0.  ]
 [0.   0.33 0.   0.33 0.   0.   0.   0.33 0.   0.  ]
 [0.   0.   0.33 0.   0.33 0.   0.   0.   0.33 0.  ]
 [0.   0.   0.   0.5  0.   0.   0.   0.   0.   0.5 ]
 [0.33 0.   0.   0.   0.   0.   0.33 0.   0.   0.  ]
 [0.   0.25 0.   0.   0.   0.25 0.   0.25 0.   0.  ]
 [0.   0.   0.25 0.   0.   0.   0.25 0.   0.25 0.  ]
 [0.   0.   0.   0.25 0.   0.   0.   0.25 0.   0.25]
 [0.   0.   0.   0.   0.33 0.   0.   0.   0.33 0.  ]]

　空間隣接行列を行ごとに基準化(行和で割る)して、空間重み行列を作成します。

　フレームごとの変数の変化を設定します。

# 変数の値を指定
x_val = 10.0

# 追加する区域番号を指定
if n%2 == 1: # 区域数が奇数の場合
    var_add_i  = [np.arange(0, N, step=2)] # 互い違いのインデックス
    var_add_i += [np.arange(1, N, step=2)] # 隙間のインデックス
    var_add_i  = np.hstack(var_add_i)
else: # 区域数が偶数の場合
    var_add_i  = [(n*i)+np.arange(i%2, n, step=2) for i in range(n)] # 互い違いのインデックス
    var_add_i += [(n*i)+np.arange((i+1)%2, n, step=2) for i in range(n)] # 隙間のインデックス
    var_add_i  = np.array(var_add_i)
print(var_add_i)

# 除去する区域番号を指定
var_del_i = np.arange(N)
print(var_del_i)

[ 0  2  4  6  8 10 12 14 16 18 20 22 24  1  3  5  7  9 11 13 15 17 19 21
 23]
[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 24]

　この例では、全ての区域が値を持たない(値が 0 の)状態を初期値として、値を持つ区域を隣接しないように増やしていき、続いて全ての区域が値を持つようにします。さらに、順番に値を消していきます。
　各区域の変数の値を x_val として指定しておきます。

　GMIの計算のヒートマップのアニメーションを作成します。

# フレーム数を取得
frame_num = len(np.hstack([var_add_i, var_del_i]))

# 色の設定用
x_max          = abs(x_val)
x_dev_max      = abs(x_val)
#x_dev_max      = abs(x_val**2)
x_dev_prod_max = abs(x_val**2)
w_max  = 1.0
z_max  = 0.15
_z_max_lt = [] # 確認用

# ラベル用の文字列を作成
fnc_label   = '$\\frac{N}{\\sum_{i=1}^N \\sum_{j=1}^N w_{ij}} '
fnc_label  += 'W \\odot (x-\\bar{x}) (x-\\bar{x})^T'
fnc_label  += '\\frac{1}{\\sum_{i=1}^N (x_i-\\bar{x})}$'
elem_label  = '$\\frac{N}{\\sum_{i=1}^N \\sum_{j=1}^N w_{ij}} '
elem_label += '\\frac{w_{ij} (x_i-\\bar{x}) (x_j-\\bar{x})}{\\sum_{i=1}^N (x_i-\\bar{x})}$'

# サイズ調整用の値を指定
d = 0.8

# 変数を初期化
x = np.zeros(N)
city_mat = np.tile(np.nan, reps=(n+2, n+2))

# グラフオブジェクトを初期化
fig, axes = plt.subplots(nrows=2, ncols=4, constrained_layout=True, 
                         figsize=((12+2+2+12)*d, (10+10)*d), width_ratios=[10, 2, 2, 10], 
                         dpi=100, facecolor='white')
fig.suptitle("Global Moran's I", fontsize=20)

# 作図処理を定義
def update(frame_i):
    
    # 前フレームのグラフを初期化
    [ax.cla() for ax in axes.flatten()]

    # 0除算の回避用
    eps = np.random.normal(loc=0.0, scale=0.001, size=1).item()

    # 変数を変更
    if frame_i < N:
        # 値を追加
        var_i     = var_add_i[frame_i]
        x[var_i]  = x_val
        x[var_i] += eps
    else:
        # 値を除去
        var_i    = var_del_i[frame_i-N]
        x[var_i] = eps

    # 格子状の区域に変形
    x_mat = x.reshape((n, n))
    city_mat[1:(n+1), 1:(n+1)] = x_mat
    
    # 偏差を計算
    x_dev = (x - np.mean(x))[:, np.newaxis]
    
    # 偏差の積を計算
    X_dev_prod = x_dev @ x_dev.T
    
    # 基準化済み重み付け偏差の積を計算
    Z = N/np.sum(W) * W*X_dev_prod / np.sum(x_dev**2)

    # 色設定用の最大値の確認用
    _z_max_lt.append(abs(Z).max())
    print('z max:', max(_z_max_lt))
    
    # 区域を作図
    ax = axes[0, 0]
    ax.pcolor(city_mat, cmap='jet', vmin=-x_max, vmax=x_max) # 変数
    city_i = 0
    for label_i in range(1, n+1):
        for label_j in range(1, n+1):
            city_i += 1
            ax.text(x=label_j+0.5, y=label_i+0.5, s=f'$i = {city_i}$', 
                    size=15, ha='center', va='center') # 区域番号
    ax.set_xticks(ticks=np.arange(1, n+1)+0.5)
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(1, n+1)+0.5)
    ax.set_yticklabels(labels='')
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('longitude')
    ax.set_ylabel('latitude')
    ax.set_title(f'$GMI = {culc_GMI(x, W):.3f}, \\bar{{x}} = {np.mean(x):.2f}$', loc='left')
    ax.grid()
    ax.set_aspect('equal', adjustable='box')
    
    # 変数を作図
    ax = axes[0, 1]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(x[:, np.newaxis], cmap='jet', vmin=-x_max, vmax=x_max) # 変数
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title('$x$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$x_i$')
    ax.set_aspect('equal', adjustable='box')
    
    # 偏差を作図
    ax = axes[0, 2]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='100%', pad=0.25)
    c = ax.pcolor(x_dev, cmap='jet', vmin=-x_dev_max, vmax=x_dev_max) # 変数
    ax.set_xticks(ticks=[0.5])
    ax.set_xticklabels(labels='')
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_ylabel('$i$')
    ax.set_title('$x - \\bar{x}$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$x_i - \\bar{x}$')
    ax.set_aspect('equal', adjustable='box')
    
    # 偏差の積を作図
    ax = axes[0, 3]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(X_dev_prod, cmap='jet', vmin=-x_dev_prod_max, vmax=x_dev_prod_max) # 変数の積
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_ylabel('$i$')
    ax.set_title('$(x-\\bar{x}) (x-\\bar{x})^T$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$(x_i-\\bar{x}) (x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')
    
    # 空間重み行列を作図
    ax = axes[1, 0]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(W, cmap='seismic', vmin=-w_max, vmax=w_max) # 重み
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_ylabel('$i$')
    ax.set_title('$W$', fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label='$(x_i-\\bar{x}) (x_j-\\bar{x})$')
    ax.set_aspect('equal', adjustable='box')
    
    # 基準化済み重み付け偏差の積を作図
    ax = axes[1, 3]
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size=f'{1/N*100}%', pad=0.25)
    c = ax.pcolor(Z, cmap='RdYlBu_r', vmin=-z_max, vmax=z_max) # 基準化済み重み付け偏差の積
    ax.set_xticks(ticks=np.arange(N)+0.5)
    ax.set_xticklabels(labels=np.arange(N)+1)
    ax.set_yticks(ticks=np.arange(N)+0.5)
    ax.set_yticklabels(labels=np.arange(N)+1)
    ax.invert_yaxis() # 軸の反転
    ax.set_xlabel('$j$')
    ax.set_ylabel('$i$')
    ax.set_title(fnc_label, fontsize=20)
    ax.grid()
    fig.colorbar(c, cax=cax, label=elem_label)
    ax.set_aspect('equal', adjustable='box')
    
    # 不要な図を非表示
    axes[1, 1].axis('off')
    axes[1, 2].axis('off')

# 動画を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=500)

# 動画を書出
ani.save(
    filename='ch2_3_gmi_adj.mp4', 
    progress_callback = lambda i, n: print(f'frame: {i} / {n}')
)

　一定の値を持つ区域の数や位置を変化させて、GMIへの影響を確認します。

　左上図は、簡易的な表現のコロプレス図(地図)です(行列を可視化した図ではありません)。上下左右の区域に隣接しています。
　中上図は、変数と偏差のベクトルです。変数の値によって平均や偏差の値に影響します。
　左下図は、隣接重み行列です。
　右上図は、隣接重み行列によって重み付けした偏差の積です。
　右下図は、基準化した重み付け偏差の積です。総和がGMIです。

　ここまでの記事では、グローバルなモランのIを確認しました。次からの記事では、ローカルなモランのIを確認します。

参考文献

作者:植杉　大
Independently published

おわりに

　重み付けして和をとる操作は加重平均の計算なわけですが、隣接重み行列の場合は、隣接していない区域に0を掛けて変数の影響を消す操作でもあるんですね。この記事の中でそれを説明すると話の流れが悪くなったので、空間ラグの記事で説明して、本での順番と変わりますがこの記事の前の記事にすることにしました。この記事を更新した時点ではまだ書き終わっていませんが。

　昨日に引き続きただの日記ですが、今日もえびちゅうのライブに行ってきましたー。今日は千秋楽でした。地元会場かつ千秋楽の公演だったので半年ほど前に申し込んだのですが、行って良かったー。ラストを大阪にしてくれてありがとうございます！
　とういうわけでまたまたえびちゅう曲を置いておくので聴いてください。

　次は半年後かなぁ。それまで頑張る。

【次の内容】