はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学時のまとめノートです。「定義編・導出編」「作図編・実装編」「可視化編」の三部構成でモデルやアルゴリズムの理解を目指します。
　本の内容から寄り道・回り道しながら進めます。本を読んだ上で補助的に読んでください。

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2.3 グローバル・モランのIの実装
参考文献
おわりに

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

　空間的自己相関(spatial autocorrelation)の指標であるグローバル・モランのI(GMI・global Moran's I)の計算を実装します。
　GMIについては「2.3：グローバル・モランのIの定義式【空間DS入門のノート】 - からっぽのしょこ」を参照してください。

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

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

疑似データの生成

　まずは、地理空間データを模した簡易的なデータ(トイデータ)を作成します。
　空間隣接行列については「【Python】2.1：空間隣接行列の作図：境界の共有【空間データサイエンス入門のノート】 - からっぽのしょこ」、空間重み行列については「【Python】2.2：空間重み行列の作図：境界の共有【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　入力データを作成します。

# 地域数を指定
N = 10

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

[1.22 0.28 0.25 0.85 3.81 0.68 0.14 0.5  3.2  0.81]
(10,)

　 $N$ 個の地域のデータ $\mathbf{x}$ を作成します。

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

　この例では、対数公示価格などの疑似データとして、対数正規乱数を使います。
　対数正規分布の乱数は、np.random.lognormal() で生成できます。サンプルサイズの引数 size に $N$ を指定します。

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

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

# 対称行列に変換
random_mat *= random_mat.T

# 空間隣接行列に変換
adj_mat = (random_mat > 0).astype(dtype=np.int8) # バイナリ値に変換
adj_mat[np.arange(N), np.arange(N)] = 0          # 対角要素を0に置換
print(adj_mat)
print(adj_mat.shape)

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

[[0 1 1 0 1 0 1 1 1 0]
 [1 0 1 0 0 0 0 1 1 0]
 [1 1 0 1 0 1 1 1 0 0]
 [0 0 1 0 1 0 1 0 1 0]
 [1 0 0 1 0 1 1 0 1 1]
 [0 0 1 0 1 0 1 1 0 1]
 [1 0 1 1 1 1 0 1 1 0]
 [1 1 1 0 0 1 1 0 0 0]
 [1 1 0 1 1 0 1 0 0 0]
 [0 0 0 0 1 1 0 0 0 0]]
(10, 10)
[[0.   0.17 0.17 0.   0.17 0.   0.17 0.17 0.17 0.  ]
 [0.25 0.   0.25 0.   0.   0.   0.   0.25 0.25 0.  ]
 [0.17 0.17 0.   0.17 0.   0.17 0.17 0.17 0.   0.  ]
 [0.   0.   0.25 0.   0.25 0.   0.25 0.   0.25 0.  ]
 [0.17 0.   0.   0.17 0.   0.17 0.17 0.   0.17 0.17]
 [0.   0.   0.2  0.   0.2  0.   0.2  0.2  0.   0.2 ]
 [0.14 0.   0.14 0.14 0.14 0.14 0.   0.14 0.14 0.  ]
 [0.2  0.2  0.2  0.   0.   0.2  0.2  0.   0.   0.  ]
 [0.2  0.2  0.   0.2  0.2  0.   0.2  0.   0.   0.  ]
 [0.   0.   0.   0.   0.5  0.5  0.   0.   0.   0.  ]]
(10, 10)

　空間重み行列 $\mathbf{W}$ を作成します。

$\displaystyle \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}$

　この例では、一様乱数を使って、空間隣接行列を作成します。
　連続一様分布の乱数は、np.random.uniform() で作成できます。確率変数の範囲の引数 low, high に0を中心とする範囲、サンプルサイズの引数 size に $N \times N$ を指定します。
　乱数の行列に転置した行列を掛けて対称行列にします。
　因子型の値 True, False は整数型の値 1, 0 として扱えることから、対称行列の乱数の正の値を 1、負の値を 0 に変換して、空間隣接行列とします。
　空間隣接行列を行ごとに正規化(規準化)して、空間重み行列を作成します。

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

スポンサードリンク

統計量の計算

　次は、グローバル・モランのIの計算を行う処理を確認します。

1次元配列の場合

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

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

# 空間重み行列を設定
W = weight_mat.copy()

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

[1.22 0.28 0.25 0.85 3.81 0.68 0.14 0.5  3.2  0.81]
(10,)
10

　入力データ $\mathbf{x}$ を要素数 $N$ の1次元配列(横ベクトル)とします。

　要素ごとの計算により、GMIを計算します。

# 標本平均を計算
#x_bar = np.sum(x) / N
x_bar = np.mean(x)
print(x_bar.round(2))

# 偏差を計算
x_dev = x - x_bar
print(x_dev.round(2))
print(x_dev.shape)

# 標本分散を計算
#x_var = np.sum(x_dev**2) / N
x_var = np.var(x, ddof=0)
print(x_var.round(2))

# GMIを計算:(外積)
GMI  = np.sum(
    W * np.outer(x_dev, x_dev)
)
GMI /= x_var
GMI /= np.sum(W)
print(GMI)

1.17
[ 0.04 -0.89 -0.92 -0.33  2.64 -0.5  -1.03 -0.67  2.03 -0.37]
(10,)
1.47
0.0037935822724658305

　GMIを、次の式で計算します。

$\displaystyle \begin{aligned} \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 } \\ &= \frac{1}{S_0} \frac{1}{s^2} \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x}) \end{aligned}$

　ここで、 $\bar{x}$ は変数 $\mathbf{x}$ の標本平均、 $s^2$ は変数 $\mathbf{x}$ の標本分散です。また、空間重み行列 $\mathbf{W}$ の総和を $S_0 = \sum_{i=1}^N \sum_{j=1}^N w_{ij}$ とおきます。

　偏差ベクトルの外積(2つのベクトルの要素の全ての組み合わせの積の行列)を np.outer() で計算し、空間重み行列の同じインデックスの要素との積を計算して、GMIを求めます。
　または、ブロードキャストにより、GMIを計算します。

# GMIを計算:(ブロードキャスト)
GMI  = np.sum(
    x_dev.reshape((N, 1)) * W * x_dev
)
GMI /= x_var
GMI /= np.sum(W)
print(GMI)

0.0037935822724658084

　空間重み行列の各行と各列に偏差ベクトルを掛けて、総和を求めます。

2次元配列の場合

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

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

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

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

[1.22 0.28 0.25 0.85 3.81 0.68 0.14 0.5  3.2  0.81]
(10, 1)
10

　入力データ $\mathbf{x}$ を $N \times 1$ の2次元配列(縦ベクトル)とします。

　行列演算により、GMIを計算します。

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

# 偏差を計算
x_dev = x - x_bar
print(x_dev.round(2).flatten())
print(x_dev.shape)

# 標本分散を計算:(行列演算)
x_var = (x_dev.T @ x_dev / N).squeeze()
print(x_var.round(2))

# GMIを計算:(行列演算)
GMI  = (x_dev.T @ W @ x_dev).squeeze()
GMI /= x_var
GMI /= np.sum(W)
print(GMI)

1.17
[ 0.04 -0.89 -0.92 -0.33  2.64 -0.5  -1.03 -0.67  2.03 -0.37]
(10, 1)
1.47
0.0037935822724657897

　GMIを、次の式で計算します。

$\displaystyle \begin{aligned} s^2 &= \frac{1}{N} \tilde{\mathbf{x}}^{\top} \tilde{\mathbf{x}} \\ \mathrm{GMI} &= \frac{1}{S_0} \frac{1}{s^2} \tilde{\mathbf{x}}^{\top} \mathbf{W} \tilde{\mathbf{x}} \end{aligned}$

　ここで、偏差を $\tilde{x}_n = x_n - \bar{x}$ 、偏差ベクトルを $\tilde{\mathbf{x}}^{\top} = (\tilde{x}_1, \cdots, \tilde{x}_N)$ で表しています。

　@ 演算子を使って、偏差ベクトルの内積(要素ごとの積の和)で標本分散、偏差ベクトルと空間重み行列の二次形式を計算して、GMIを求めます。

　入力データ x が1次元配列の場合もこのコードで処理できます。

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

スポンサードリンク

統計量の実装

　続いて、グローバル・モランのIの計算を行う関数を作成します。

実装

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

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

    # 分散を計算
    x_var = np.var(x, ddof=ddof)
    #x_var = (x_dev.T @ x_dev / (N - ddof)).squeeze()
    
    # GMIを計算
    GMI  = (x_dev.T @ W @ x_dev).squeeze()
    GMI /= np.sum(W) * x_var
    return GMI

　「2次元配列の場合」の処理を関数として実装します。
　入力データ x が1次元・2次元どちらの場合も処理できます。入力データ x の形状に関わらず、出力データ GMI はスカラなります。

計算

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

# データを整形
x = x.flatten()
print(x.shape)

# GMIを計算
GMI = calc_GMI(x, W)
print(GMI)
print(GMI.shape)

(10,)
0.0037935822724657893
()

# データを整形
x = x.reshape((N, 1))
print(x.shape)

# GMIを計算
GMI = calc_GMI(x, W)
print(GMI)
print(GMI.shape)

(10, 1)
0.0037935822724657893
()

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

　この記事では、GMIをプログラムで確認しました。次の記事では、優位性検定をプログラムで確認します。

参考文献

Pythonで学ぶ空間データサイエンス入門: 地域の特徴を発見する方法

作者:植杉　大
Independently published

Amazon

おわりに

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

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

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

2026.07.30：加筆修正しました。

　他の記事の修正に合わせて解説のレベル感や表現を調整しました。

　その後も何度もってほどではないですが何度か単独やツーマン、フェスでえびちゅうのライブに行けてます。やっぱりライブは良いですねぇ。ところで登校って設定のことはすっかり忘れてました。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】2.3：グローバル・モランのIの実装【空間DS入門のノート】