はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

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

2.4 ローカル・モランのIの実装

　空間的自己相関(spatial autocorrelation)の指標であるローカル・モランのI(LMI・local Moran's I)の計算を実装します。
　LMIについては「2.4：ローカル・モランのIの定義式【空間DS入門のノート】 - からっぽのしょこ」、空間ラグ(SL・spatial lag)については「【Python】2.3：空間ラグの実装【空間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)

[0.9  2.64 0.98 1.15 1.14 1.34 1.03 2.14 2.95 2.75]
(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 1 0 0 0 0 0 1]
 [1 0 1 0 0 0 0 0 1 0]
 [1 1 0 1 0 0 1 0 0 1]
 [1 0 1 0 1 0 0 0 0 0]
 [0 0 0 1 0 0 1 1 1 1]
 [0 0 0 0 0 0 0 0 1 1]
 [0 0 1 0 1 0 0 0 1 0]
 [0 0 0 0 1 0 0 0 0 0]
 [0 1 0 0 1 1 1 0 0 1]
 [1 0 1 0 1 1 0 0 1 0]]
(10, 10)
[[0.   0.25 0.25 0.25 0.   0.   0.   0.   0.   0.25]
 [0.33 0.   0.33 0.   0.   0.   0.   0.   0.33 0.  ]
 [0.2  0.2  0.   0.2  0.   0.   0.2  0.   0.   0.2 ]
 [0.33 0.   0.33 0.   0.33 0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.2  0.   0.   0.2  0.2  0.2  0.2 ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.5  0.5 ]
 [0.   0.   0.33 0.   0.33 0.   0.   0.   0.33 0.  ]
 [0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  ]
 [0.   0.2  0.   0.   0.2  0.2  0.2  0.   0.   0.2 ]
 [0.2  0.   0.2  0.   0.2  0.2  0.   0.   0.2  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 に変換して、空間隣接行列とします。
　空間隣接行列を行ごとに正規化(規準化)して、空間重み行列を作成します。

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

スポンサードリンク

統計量の計算

　次は、ローカル・モランのIの計算における処理を確認します。

1次元配列の場合

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

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

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

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

[0.9  2.64 0.98 1.15 1.14 1.34 1.03 2.14 2.95 2.75]
(10,)
10

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

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

# 標本平均を計算
#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))

# LMIを計算:(ループ処理)
SL  = np.tile(np.nan, reps=N) # 受け皿を初期化
LMI = np.tile(np.nan, reps=N) # 受け皿を初期化
for i in range(N):
    SL[i]  = np.sum(W[i] * x_dev)     # 値を格納
    LMI[i] = x_dev[i] / x_var * SL[i] # 値を格納
print(LMI.round(2))
print(LMI.shape)

1.7
[-0.8   0.94 -0.72 -0.55 -0.56 -0.37 -0.67  0.43  1.25  1.05]
(10,)
0.61
[-0.24 -0.14  0.01  0.63 -0.28 -0.69  0.01 -0.4   0.16 -0.41]
(10,)

　各地域のLMIは、次の式で計算できます。

$\displaystyle \begin{aligned} \mathrm{LMI}_i &= \frac{ (x_i - \bar{x}) \sum_{j=1}^N w_{ij} (x_j - \bar{x}) }{ \frac{1}{N} \sum_{j=1}^N (x_j - \bar{x})^2 } \\ &= \frac{x_i - \bar{x}}{s^2} \mathrm{SL}_i \end{aligned}$

　ここで、 $\bar{x}$ は変数 $\mathbf{x}$ の標本平均、 $s^2$ は変数 $\mathbf{x}$ の標本分散です。
　標本分散 $s^2 = \frac{1}{N} \sum_{n=1}^N (x_n - \bar{x})^2$ を用いる場合の定義式です。不偏分散 $\hat{\sigma}^2 = \frac{1}{N-1} \sum_{n=1}^N (x_n - \bar{x})^2$ を用いて定義する場合もあります。

　for 文を使って、順番に地域ごとのLMIを求めて配列に格納します。
　または、ブロードキャストにより、LMIを計算します。

# 空間ラグを計算:(ブロードキャスト)
SL = np.sum(W * x_dev, axis=1)
print(SL.round(2))
print(SL.shape)

# LMIを計算:(ブロードキャスト)
LMI  = x_dev / x_var * SL
print(LMI.round(2))
print(LMI.shape)

[ 0.18 -0.09 -0.01 -0.69  0.3   1.15 -0.01 -0.56  0.08 -0.24]
(10,)
[-0.24 -0.14  0.01  0.63 -0.28 -0.69  0.01 -0.4   0.16 -0.41]
(10,)

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

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)

[0.9  2.64 0.98 1.15 1.14 1.34 1.03 2.14 2.95 2.75]
(10, 1)
10

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

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

# 標本平均を計算
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))

# 空間ラグを計算:(行列演算)
SL = W @ x_dev
print(SL.round(2).flatten())
print(SL.shape)

# LMIを計算:(行列演算)
LMI  = x_dev / x_var * SL
print(LMI.round(2).flatten())
print(LMI.shape)

1.7
[-0.8   0.94 -0.72 -0.55 -0.56 -0.37 -0.67  0.43  1.25  1.05]
(10, 1)
0.61
[ 0.18 -0.09 -0.01 -0.69  0.3   1.15 -0.01 -0.56  0.08 -0.24]
(10, 1)
[-0.24 -0.14  0.01  0.63 -0.28 -0.69  0.01 -0.4   0.16 -0.41]
(10, 1)

　 $N$ 個の地域のLMIを、次の式で計算します。

$\displaystyle \begin{aligned} s^2 &= \frac{1}{N} \tilde{\mathbf{x}}^{\top} \tilde{\mathbf{x}} \\ \begin{pmatrix} \mathrm{SL}_1 \\ \mathrm{SL}_2 \\ \vdots \\ \mathrm{SL}_N \end{pmatrix} &= \mathbf{W} \tilde{\mathbf{x}} \\ \begin{pmatrix} \mathrm{LMI}_1 \\ \mathrm{LMI}_2 \\ \vdots \\ \mathrm{LMI}_N \end{pmatrix} &= \frac{1}{s^2} \tilde{\mathbf{x}} \odot \mathbf{W} \tilde{\mathbf{x}} \end{aligned}$

　ここで、 $\odot$ はアダマール積(要素ごとの積)の演算子です。また、偏差を $\tilde{x}_n = x_n - \bar{x}$ 、偏差ベクトルを $\tilde{\mathbf{x}}^{\top} = (\tilde{x}_1, \cdots, \tilde{x}_N)$ で表しています。

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

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

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

スポンサードリンク

統計量の実装

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

実装

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

# LMIの計算を実装
def calc_LMI(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()

    # SLを計算
    SL = W @ x_dev

    # LMIを計算
    LMI  = x_dev / x_var * SL
    return LMI

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

計算

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

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

# LMIを計算
LMI = calc_LMI(x, W)
print(LMI.round(2).flatten())
print(LMI.shape)

(10,)
[-0.24 -0.14  0.01  0.63 -0.28 -0.69  0.01 -0.4   0.16 -0.41]
(10,)

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

# LMIを計算
LMI = calc_LMI(x, W)
print(LMI.round(2).flatten())
print(LMI.shape)

(10, 1)
[-0.24 -0.14  0.01  0.63 -0.28 -0.69  0.01 -0.4   0.16 -0.41]
(10, 1)

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

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

参考文献

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

作者:植杉　大
Independently published

Amazon

おわりに

　抽象的な概念を視覚的に理解できるようにしようと書き始めたシリーズですが、結局いつものように数式から数式を見る記事、数式をプログラミングして再現する記事、プログラムを図にして数式との対応を見る記事の3タイプの構成になりそうです。というわけで次の記事に続きます。

2026.07.30：加筆修正しました。

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

　7月30日は、モーニング娘。の6期メンバーのCDデビューです。

　ちゃゆ。

【次の内容】

からっぽのしょこ

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

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