はじめに

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

　この記事では、空間ラグについて、プログラムを使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2.3 空間ラグの実装
参考文献
おわりに

2.3 空間ラグの実装

　空間ラグ(SL・spatial lag)の計算処理を実装します。
　空間重み行列(spatial weight matrix)については「【Python】2.2：空間重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

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

空間ラグの数式表現

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

　 $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$ 番目の区域のSLは、次の式で定義されます。

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

　SLの計算に関して、 $\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}$

とすると、 $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} = \begin{pmatrix} \mathrm{SL}_1 \\ \mathrm{SL}_2 \\ \vdots \\ \mathrm{SL}_N \end{pmatrix}$

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

[0.19 1.28 0.84 0.5  0.62 1.39 1.71 1.05 0.89 1.06]
(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.98 -0.42 -0.02  0.47 -0.23 -0.31 -0.58 -0.16  0.01  0.29]
 [-0.42  0.46  0.18  0.01 -0.11 -0.03  0.35 -0.04  0.12  0.4 ]
 [-0.02  0.18  0.3  -0.05  0.24  0.86 -0.38 -0.19  0.49 -0.19]
 [ 0.47  0.01 -0.05  0.38  0.02 -0.36 -0.17 -0.03  0.05 -0.24]
 [-0.23 -0.11  0.24  0.02  0.34  0.68 -0.15  0.06  0.03  0.03]
 [-0.31 -0.03  0.86 -0.36  0.68  0.04 -0.24 -0.01 -0.05 -0.55]
 [-0.58  0.35 -0.38 -0.17 -0.15 -0.24  0.23 -0.04  0.04 -0.75]
 [-0.16 -0.04 -0.19 -0.03  0.06 -0.01 -0.04  0.04 -0.47 -0.14]
 [ 0.01  0.12  0.49  0.05  0.03 -0.05  0.04 -0.47  0.    0.06]
 [ 0.29  0.4  -0.19 -0.24  0.03 -0.55 -0.75 -0.14  0.06  0.11]]
[[0 0 0 1 0 0 0 0 1 1]
 [0 0 1 1 0 0 1 0 1 1]
 [0 1 0 0 1 1 0 0 1 0]
 [1 1 0 0 1 0 0 0 1 0]
 [0 0 1 1 0 1 0 1 1 1]
 [0 0 1 0 1 0 0 0 0 0]
 [0 1 0 0 0 0 0 0 1 0]
 [0 0 0 0 1 0 0 0 0 0]
 [1 1 1 1 1 0 1 0 0 1]
 [1 1 0 0 1 0 0 0 1 0]]
[[0.   0.   0.   0.33 0.   0.   0.   0.   0.33 0.33]
 [0.   0.   0.2  0.2  0.   0.   0.2  0.   0.2  0.2 ]
 [0.   0.25 0.   0.   0.25 0.25 0.   0.   0.25 0.  ]
 [0.25 0.25 0.   0.   0.25 0.   0.   0.   0.25 0.  ]
 [0.   0.   0.17 0.17 0.   0.17 0.   0.17 0.17 0.17]
 [0.   0.   0.5  0.   0.5  0.   0.   0.   0.   0.  ]
 [0.   0.5  0.   0.   0.   0.   0.   0.   0.5  0.  ]
 [0.   0.   0.   0.   1.   0.   0.   0.   0.   0.  ]
 [0.14 0.14 0.14 0.14 0.14 0.   0.14 0.   0.   0.14]
 [0.25 0.25 0.   0.   0.25 0.   0.   0.   0.25 0.  ]]
(10, 10)

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

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

空間ラグのIの計算

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

1次元配列で扱う場合

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

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

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

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

[0.19 1.28 0.84 0.5  0.62 1.39 1.71 1.05 0.89 1.06]
(10,)

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

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

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

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

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

0.95
[-0.77  0.33 -0.11 -0.45 -0.33]
[-0.14  0.05  0.09 -0.21  0.   -0.22  0.13 -0.33 -0.07 -0.21]
(10,)

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

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

2次元配列で扱う場合

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

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

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

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

[[0.19]
 [1.28]
 [0.84]
 [0.5 ]
 [0.62]
 [1.39]
 [1.71]
 [1.05]
 [0.89]
 [1.06]]
(10, 1)

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

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

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

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

# SLを計算
SL = W @ x_dev
print(SL.round(2))
print(SL.shape)
print(SL.flatten().round(2)) # 1次元配列に変換
print(SL.flatten().shape)

0.95
[[-0.77]
 [ 0.33]
 [-0.11]
 [-0.45]
 [-0.33]]
[[-0.14]
 [ 0.05]
 [ 0.09]
 [-0.21]
 [ 0.  ]
 [-0.22]
 [ 0.13]
 [-0.33]
 [-0.07]
 [-0.21]]
(10, 1)
[-0.14  0.05  0.09 -0.21  0.   -0.22  0.13 -0.33 -0.07 -0.21]
(10,)

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

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

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

空間ラグの実装

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

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

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

    # SLを計算
    SL = W @ x_dev
    return SL.flatten()

　入力変数 x が (N, 1) の2次元配列(縦ベクトル)の場合、出力変数 SL が (N, 1) の2次元配列になるので、.flatten() で1次元配列として出力します。
　x が (N,) の1次元配列の場合もそのまま処理できます。

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

# SLを計算
SL = culc_SL(x, W)
print(SL.round(2))

[-0.14  0.05  0.09 -0.21  0.   -0.22  0.13 -0.33 -0.07 -0.21]

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

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

参考文献

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

作者:植杉　大
Independently published

Amazon

おわりに

　空間ラグ自体はGMIに出てきませんが、簡単例として先に頭に入れておいた方が諸々をイメージしやすいと思うので、こちらを先に扱っておくことにします。

【次の内容】

からっぽのしょこ

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

【Python】2.3：空間ラグの実装【空間データサイエンス入門のノート】