はじめに

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

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

【前の内容】

【他の内容】

【今回の内容】

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

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

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

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

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

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

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

　 $i$ 番目の区域のLMIは、次の式で定義されます。

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

　不偏分散を用いる場合の定義式です。

　LMIの計算に関して、 $\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_{j=1}^N (x_i - \bar{x})^2$

　また、 $N$ 個の区域に関する重み付けした偏差の和は、次の行列とベクトルの積で計算できます。

$\displaystyle \mathbf{W} \tilde{\mathbf{x}} = \begin{pmatrix} \sum_{j=1}^N w_{1j} (x_j - \bar{x}) \\ \sum_{j=1}^N w_{2j} (x_j - \bar{x}) \\ \vdots \\ \sum_{j=1}^N w_{Nj} (x_j - \bar{x}) \end{pmatrix}$

　 $N$ 個の区域に関する偏差と重み付け偏差の和の積は、ベクトルの要素ごとの積で計算できます。

$\displaystyle \tilde{\mathbf{x}} \odot \mathbf{W} \tilde{\mathbf{x}} = \begin{pmatrix} (x_1 - \bar{x}) \sum_{j=1}^N w_{1j} (x_j - \bar{x}) \\ (x_2 - \bar{x}) \sum_{j=1}^N w_{2j} (x_j - \bar{x}) \\ \vdots \\ (x_N - \bar{x}) \sum_{j=1}^N w_{Nj} (x_j - \bar{x}) \end{pmatrix}$

　 $\odot$ はアダマール積の演算子であり、対応する要素ごとの積を表します。

　式変形については「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)

[1.   3.67 2.13 0.65 1.63 0.44 1.34 5.13 0.29 0.59]
(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.16 -0.2  -0.37  0.33  0.26  0.01 -0.44  0.11 -0.48  0.12]
 [-0.2   0.68 -0.74 -0.02 -0.02 -0.01  0.05 -0.27 -0.3  -0.29]
 [-0.37 -0.74  0.4  -0.13  0.04 -0.09 -0.2  -0.35  0.01  0.32]
 [ 0.33 -0.02 -0.13  0.    0.34 -0.19  0.69  0.1  -0.62  0.16]
 [ 0.26 -0.02  0.04  0.34  0.63  0.08  0.54  0.17  0.06  0.09]
 [ 0.01 -0.01 -0.09 -0.19  0.08  0.53 -0.13  0.1  -0.1   0.02]
 [-0.44  0.05 -0.2   0.69  0.54 -0.13  0.   -0.39 -0.07  0.08]
 [ 0.11 -0.27 -0.35  0.1   0.17  0.1  -0.39  0.04 -0.53 -0.  ]
 [-0.48 -0.3   0.01 -0.62  0.06 -0.1  -0.07 -0.53  0.01  0.5 ]
 [ 0.12 -0.29  0.32  0.16  0.09  0.02  0.08 -0.    0.5   0.36]]
[[0 0 0 1 1 1 0 1 0 1]
 [0 0 0 0 0 0 1 0 0 0]
 [0 0 0 0 1 0 0 0 1 1]
 [1 0 0 0 1 0 1 1 0 1]
 [1 0 1 1 0 1 1 1 1 1]
 [1 0 0 0 1 0 0 1 0 1]
 [0 1 0 1 1 0 0 0 0 1]
 [1 0 0 1 1 1 0 0 0 0]
 [0 0 1 0 1 0 0 0 0 1]
 [1 0 1 1 1 1 1 0 1 0]]
[[0.   0.   0.   0.2  0.2  0.2  0.   0.2  0.   0.2 ]
 [0.   0.   0.   0.   0.   0.   1.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.33 0.   0.   0.   0.33 0.33]
 [0.2  0.   0.   0.   0.2  0.   0.2  0.2  0.   0.2 ]
 [0.12 0.   0.12 0.12 0.   0.12 0.12 0.12 0.12 0.12]
 [0.25 0.   0.   0.   0.25 0.   0.   0.25 0.   0.25]
 [0.   0.25 0.   0.25 0.25 0.   0.   0.   0.   0.25]
 [0.25 0.   0.   0.25 0.25 0.25 0.   0.   0.   0.  ]
 [0.   0.   0.33 0.   0.33 0.   0.   0.   0.   0.33]
 [0.14 0.   0.14 0.14 0.14 0.14 0.14 0.   0.14 0.  ]]
(10, 10)

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

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

ローカル・モランのIの計算

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

1次元配列で扱う場合

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

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

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

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

[1.   3.67 2.13 0.65 1.63 0.44 1.34 5.13 0.29 0.59]
(10,)

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

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

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

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

# LMIを計算
LMI = np.tile(np.nan, reps=N)
for i in range(N):
    LMI_i  = x_dev[i] * np.sum(W[i] * x_dev) * (N-1) / np.sum(x_dev**2)
    LMI[i] = LMI_i
print(LMI.round(2))
print(LMI.shape)

1.69
[-0.69  1.99  0.44 -1.04 -0.06]
[-0.   -0.28 -0.15 -0.1   0.01 -0.2   0.01 -1.05  0.13  0.27]
(10,)

　値を持たない( N 個の np.nan を格納した)1次元配列 LMI を作成しておき、各区域のLMIを順番に求めて格納していきます。
　重み付けした偏差の和の計算は、i 行目の隣接重み行列に偏差ベクトルを掛けて総和を求めます。さらに、i 個目の偏差と掛けます。
　偏差の2乗和の計算は、偏差ベクトルの各要素を2乗してから総和を求めます。

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

2次元配列で扱う場合

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

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

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

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

[[1.  ]
 [3.67]
 [2.13]
 [0.65]
 [1.63]
 [0.44]
 [1.34]
 [5.13]
 [0.29]
 [0.59]]
(10, 1)

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

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

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

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

# LMIを計算
LMI = (x_dev * (W @ x_dev)) * (N-1) / (x_dev.T @ x_dev)
print(LMI.round(2))
print(LMI.shape)
print(LMI.flatten().round(2)) # 1次元配列に変換
print(LMI.flatten().shape)

1.69
[[-0.69]
 [ 1.99]
 [ 0.44]
 [-1.04]
 [-0.06]]
[[-0.  ]
 [-0.28]
 [-0.15]
 [-0.1 ]
 [ 0.01]
 [-0.2 ]
 [ 0.01]
 [-1.05]
 [ 0.13]
 [ 0.27]]
(10, 1)
[-0.   -0.28 -0.15 -0.1   0.01 -0.2   0.01 -1.05  0.13  0.27]
(10,)

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

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

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

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

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

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

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

    # LMIを計算
    LMI = (x_dev * (W @ x_dev)) * (N-1) / (x_dev.T @ x_dev)
    return LMI.flatten()

　入力変数 x が (N, 1) の2次元配列(縦ベクトル)の場合、.T で (1, N) の2次元配列(横ベクトル)に転置して、内積を計算します。x が縦ベクトルの場合、出力変数 LMI が (N, 1) の2次元配列になるので、.flatten() で1次元配列として出力します。
　x が (N,) の1次元配列の場合もそのまま処理できます。

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

# LMIを計算
LMI = culc_LMI(x, W)
print(LMI.round(2))

[-0.   -0.28 -0.15 -0.1   0.01 -0.2   0.01 -1.05  0.13  0.27]

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

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

参考文献

作者:植杉　大
Independently published

おわりに

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

【次の内容】