はじめに

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

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

【前の内容】

【他の内容】

【今回の内容】

2.3 空間的自己相関(Global Moran's I)の実装

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

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

# ライブラリを読込
import numpy as np

グローバル・モランの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 }$

　GMIの計算に関して、 $\mathbf{x}$ の標本平均を $\bar{x}$ として、偏差(各データ $x_i$ と平均 $\bar{x}$ の差)のベクトルを

$\displaystyle \tilde{\mathbf{x}} = \begin{pmatrix} x_1 - \bar{x} \\ x_2 - \bar{x} \\ \vdots \\ x_N - \bar{x} \end{pmatrix}$

とすると、偏差の2乗和は、次の内積で計算できます。

$\displaystyle \tilde{\mathbf{x}}^{\top} \tilde{\mathbf{x}} = \sum_{i=1}^N (x_i - \bar{x})^2$

　また、重み付けした偏差の積和は、次の二次形式で計算できます。

$\displaystyle \tilde{\mathbf{x}}^{\top} \mathbf{W} \tilde{\mathbf{x}} = \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x})$

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

　以上の計算を実装します。

疑似データの生成

　実装の前に、空間データを模したデータを人工的に作成します。

　乱数を用いて、入力データを作成します。

# 区域数を指定
N = 10

# 空間データを生成
data_vec = np.random.lognormal(mean=0.0, sigma=1.0, size=N)
print(data_vec.round(2))
print(data_vec.shape)

[2.2  0.18 2.77 0.92 0.49 1.52 1.06 1.54 0.17 1.3 ]
(10,)

　この例では、対数正規分布の乱数を np.random.lognormal() で作成して、対数公示価格などのデータとして使います。

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

# 乱数を生成
random_mat = np.random.uniform(low=-1.0, high=1.0, size=(N, N))

# 対称行列に変換
random_mat *= random_mat.T
print(random_mat.round(2))

# 空間隣接行列に変換
adj_mat = (random_mat >= 0).astype('int64') # 0, 1の2値に変換
adj_mat[np.arange(N), np.arange(N)] = 0 # 対角要素を0に置換
print(adj_mat.round(2))

# 空間重み行列に変換
weight_mat = adj_mat.astype('float64')
weight_mat /= weight_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.38  0.14 -0.16 -0.57  0.34 -0.41  0.44  0.1   0.17  0.02]
 [ 0.14  0.41 -0.01  0.07  0.    0.22  0.22  0.03  0.11  0.28]
 [-0.16 -0.01  0.31  0.19 -0.05  0.01  0.17  0.6  -0.05  0.55]
 [-0.57  0.07  0.19  0.15  0.26  0.26  0.35 -0.02 -0.01 -0.06]
 [ 0.34  0.   -0.05  0.26  0.    0.65 -0.35 -0.05 -0.36  0.04]
 [-0.41  0.22  0.01  0.26  0.65  0.05 -0.24  0.66  0.07  0.03]
 [ 0.44  0.22  0.17  0.35 -0.35 -0.24  0.03  0.22  0.   -0.43]
 [ 0.1   0.03  0.6  -0.02 -0.05  0.66  0.22  0.16 -0.06 -0.21]
 [ 0.17  0.11 -0.05 -0.01 -0.36  0.07  0.   -0.06  0.1   0.09]
 [ 0.02  0.28  0.55 -0.06  0.04  0.03 -0.43 -0.21  0.09  0.4 ]]
[[0 1 0 0 1 0 1 1 1 1]
 [1 0 0 1 1 1 1 1 1 1]
 [0 0 0 1 0 1 1 1 0 1]
 [0 1 1 0 1 1 1 0 0 0]
 [1 1 0 1 0 1 0 0 0 1]
 [0 1 1 1 1 0 0 1 1 1]
 [1 1 1 1 0 0 0 1 1 0]
 [1 1 1 0 0 1 1 0 0 0]
 [1 1 0 0 0 1 1 0 0 1]
 [1 1 1 0 1 1 0 0 1 0]]
[[0.   0.17 0.   0.   0.17 0.   0.17 0.17 0.17 0.17]
 [0.12 0.   0.   0.12 0.12 0.12 0.12 0.12 0.12 0.12]
 [0.   0.   0.   0.2  0.   0.2  0.2  0.2  0.   0.2 ]
 [0.   0.2  0.2  0.   0.2  0.2  0.2  0.   0.   0.  ]
 [0.2  0.2  0.   0.2  0.   0.2  0.   0.   0.   0.2 ]
 [0.   0.14 0.14 0.14 0.14 0.   0.   0.14 0.14 0.14]
 [0.17 0.17 0.17 0.17 0.   0.   0.   0.17 0.17 0.  ]
 [0.2  0.2  0.2  0.   0.   0.2  0.2  0.   0.   0.  ]
 [0.2  0.2  0.   0.   0.   0.2  0.2  0.   0.   0.2 ]
 [0.17 0.17 0.17 0.   0.17 0.17 0.   0.   0.17 0.  ]]
(10, 10)

　形状が (N, N) の負の値を含む一様分布の乱数を np.random.uniform() で作成し、転置して掛けて、対称行列を作成します。
　因子型の値を整数型の値に変換すると True は 1、False は 0 になるのを利用して、正の値を 1、負の値 0 に変換し、空間隣接行列を作成します。
　空間隣接行列を行ごとに基準化(行和で割る)して、空間重み行列を作成します。

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

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

　次は、GMIの計算における処理を確認します。

1次元配列で扱う場合

　入力データを1次元配列として扱います。

# データを設定
x = data_vec.copy()
print(x.round(2))
print(x.shape)

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

# データ数を取得
N = len(x)

[2.2  0.18 2.77 0.92 0.49 1.52 1.06 1.54 0.17 1.3 ]
(10,)

　「疑似データの生成」で作成したデータをそのまま使います。

　要素ごとの計算によってGMIを計算します。

# 平均を計算
x_bar = np.mean(x)
print(x_bar.round(2))

# 偏差を計算
x_dev = x - x_bar
print(x_dev[:5].round(2))

# GMIを計算
GMI = N / np.sum(W) * np.sum(x_dev.reshape((N, 1)) * W * x_dev) / np.sum(x_dev**2)
print(GMI)

1.21
[ 0.98 -1.04  1.56 -0.29 -0.73]
-0.04132910619325073

　重み付けした偏差の積和の計算は、隣接重み行列の各行と列にそれぞれ偏差ベクトルを掛けて総和を求めます。
　偏差の2乗和の計算は、偏差ベクトルの各要素を2乗して総和を求めます。

　定義式に近い操作で処理しました。

2次元配列で扱う場合

　入力データを2次元配列として扱います。

# データを設定
x = data_vec.reshape((N, 1))
print(x.round(2))
print(x.shape)

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

# データ数を取得
N = len(x)

[[2.2 ]
 [0.18]
 [2.77]
 [0.92]
 [0.49]
 [1.52]
 [1.06]
 [1.54]
 [0.17]
 [1.3 ]]
(10, 1)

　(N, 1) の2次元配列(縦ベクトル)に変換します。

　行列計算によってGMIを計算します。

# 平均を計算
x_bar = np.mean(x)
print(x_bar.round(2))

# 偏差を計算
x_dev = x - x_bar
print(x_dev[:5].round(2))

# GMIを計算
GMI = N / np.sum(W) * (x_dev.T @ W @ x_dev) / (x_dev.T @ x_dev)
print(GMI)
print(GMI.item()) # スカラに変換

1.21
[[ 0.98]
 [-1.04]
 [ 1.56]
 [-0.29]
 [-0.73]]
[[-0.04132911]]
-0.04132910619325073

　重み付けした偏差の積和の計算は、@ 演算子を使って、偏差ベクトルと隣接重み行列の行列積を求めます。
　偏差の2乗和の計算は、@ 演算子を使って、偏差ベクトルの内積を求めます。

　ベクトルや行列の演算によって積と和を一度に処理しました。1次元配列の場合もこのコードで処理できます。

　以上の処理を実装します。

グローバル・モランの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()

　入力変数 x が (N, 1) の2次元配列(縦ベクトル)の場合、.T で (1, N) の2次元配列(横ベクトル)に転置して、内積や二次形式を計算します。x が縦ベクトルの場合、出力変数 GMI が (1, 1) の2次元配列になるので、.item() で値を取り出しスカラとして出力します。
　x が (N,) の1次元配列の場合もそのまま処理できます。

　実装した関数を使ってGMIの計算をします。

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

-0.04132910619325073

　以上で、GMIの計算を実装できました。

　この記事では、GMIをプログラムで確認しました。次の記事では、図で確認します。

参考文献

作者:植杉　大
Independently published

おわりに

　空間隣接行列や空間重み行列と同様に可視化の記事を1つ書くつもりだったのですが、構成の都合もあって実装(処理の確認)の記事を独立して書きました。今回のシリーズに限らずこれまでのシリーズも含めて今後は、これくらいの粒度感の方が書きやすいし読みやすいし良いのかもしれません。重複する内容が増えると思うのですが、それもまた勉強としては良いことなのかもしれません。

　ただの日記ですが、この記事の投稿日にえびちゅうこと私立恵比寿中学の単独ライブに初めて行ってきましたー。ちなみにグループ名が中学なこともありライブに行くことを登校と言ったり言わなかったりします。
　とういうわけで最後にえびちゅう曲を置いておくので聴いてください。

　色々あって(あるいは無くて)ライブに行くのは1年半ぶりだったんですが、やっぱりライブは良いなぁ何度でも行きたい。

【次の内容】