はじめに

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

　この記事では、GWRの重み行列について、プログラムを使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
4.3.2-3 GWRモデルの重み行列の実装
参考文献
おわりに

4.3.2-3 GWRモデルの重み行列の実装

　地理的加重回帰モデル(GWRモデル・geographically weighted regression model)における重み行列(weight matrix)の作成処理を実装します。
　重み関数については「【Python】4.3.2：GWRモデルの重み関数の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」、重み行列については「【Python】4.3.3：GWRモデルの重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

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

重み行列の数式表現

　まずは、GWRモデルの重み行列の定義を数式で確認します。

　 $N$ 個の地点に関して、基準地点 $i$ の重み行列を $N \times N$ の行列 $\mathbf{W}_i$ とします。

$\displaystyle \mathbf{W}_i = \begin{pmatrix} w_{i11} & w_{i12} & \cdots & w_{i1N} \\ w_{i21} & w_{i22} & \cdots & w_{i2N} \\ \vdots & \vdots & \ddots & \vdots \\ w_{iN1} & w_{iN2} & \cdots & w_{iNN} \end{pmatrix} \\ = \begin{pmatrix} w_{i11} & 0 & \cdots & 0 \\ 0 & w_{i22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & w_{iNN} \end{pmatrix}$

　異なる地点 $j, l$ の組み合わせに対応する要素(非対角要素)は $w_{ijl} = 0$ です。よって、重み行列は対角行列です。
　各地点(同じ地点の組み合わせ) $j$ に対応する要素(対角要素)は $w_{ijj}$ で表され、重み関数(カーネル)に従い、基準の地点 $i$ との距離に応じて値が決まります。基準の地点 $i$ に対応する要素は $w_{iii} = 1$ です。

　この重み行列をNumPy配列として作成します。

疑似データの生成

　重み行列の作成の前に、距離データを人工的に作成します。

　地点数を指定して、距離データを作成します。

# 地点数を指定
N = 10

# 距離を生成
dist_vec = np.random.uniform(low=0, high=100, size=N) # 範囲を指定
print(dist_vec.round(1))
print(dist_vec.shape)

[ 3.6 64.3 27.8 86.2 37.1  3.4 87.8 17.4  6.6 65.9]
(10,)

　地点数(データ数) $N$ を指定して、 $N$ 個の地点の距離 $d_{i1}, \cdots, d_{iN}$ を作成します。
　この例では、連続一様分布に従う乱数を使います。random モジュールの uniform() のサンプルサイズの引数 size に地点数 $N$ を指定します。最小値の引数 low と最大値の引数 high で乱数の範囲を指定できます。

　(緯度経度データを生成して、平面直角座標系に変換し、距離を計算したかったのですが、思っていたよりも面倒だったので諦めました。)

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

重み行列の計算

　次は、GWRモデルで利用する重み関数の定義式を確認して、重み行列の作成処理を確認します。

Moving Window型

　moving window型の重み関数(moving windowカーネル)による重み行列を作成します。

　moving windowカーネルは、次の式で定義されます。

$\displaystyle w_{ijj} = \begin{cases} 1 &\quad (d_{ij} \lt b) \\ 0 &\quad (d_{ij} \geq b) \end{cases}$

　 $d_{ij}$ は地点 $i, j$ の距離、 $b$ はバンド幅を表します。それぞれ距離なので、0以上の値 $d_{ij} \geq 0, b \geq 0$ を満たす必要があります。
　Moving Window型の重み関数は、0か1の2値 $w_{ijj} \in {0, 1}$ をとります。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定
b = 40.0

# 対角要素を計算:(moving window型)
w_diag = np.ones_like(dist_vec)
print(w_diag.round(2))

# バンド幅外の要素を置換
w_diag[dist_vec >= b] = 0.0
print(w_diag.round(2))
print(w_diag.shape)

# 重み行列を作成
W = np.diag(w_diag)
print(W.round(2))
print(W.shape)

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 0. 1. 0. 1. 1. 0. 1. 1. 0.]
(10,)
[[1. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
(10, 10)

　バンド幅 $b$ を指定して、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。全ての対角要素を $w_{ijj} = 1$ にしておき、 $d_{ij} \geq b$ の要素を $w_{ijj} = 0$ に置き換えます。非対角要素は $w_{ijl} = 0$ です。

指数型

　指数型の重み関数(指数カーネル・exponential kernel)による重み行列を作成します。

　指数カーネルは、次の式で定義されます。

$\displaystyle w_{ijj} = \exp \left( - \frac{d_{ij}}{b} \right)$

　 $d_{ij}$ は地点 $i, j$ の距離、 $b$ はバンド幅を表します。それぞれ距離なので、0以上の値 $d_{ij} \geq 0$ 、また $b$ に関しては0除算防止のため、非負の値 $b \gt 0$ を満たす必要があります。
　指数型の重み関数は、0を除く0から1の値 $0 \lt w_{ijj} \leq 1$ をとります。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定
b = 40.0

# 対角要素を計算:(指数型)
w_diag = np.exp(- dist_vec / b)
print(w_diag.round(2))
print(w_diag.shape)

# 重み行列を作成
W = np.diag(w_diag)
print(W.round(2))
print(W.shape)

[0.91 0.2  0.5  0.12 0.4  0.92 0.11 0.65 0.85 0.19]
(10,)
[[0.91 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.2  0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.5  0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.12 0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.4  0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.92 0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.11 0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.65 0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.85 0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.19]]
(10, 10)

　バンド幅 $b$ を指定して、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。対角要素は $w_{ijj} = \exp(- \frac{d_{ij}}{b})$ を計算します。非対角要素は $w_{ijl} = 0$ です。

ガウス型

　ガウス型の重み関数(ガウスカーネル・Gaussian kernel)による重み行列を作成します。

　ガウスカーネルは、次の式で定義されます。

$\displaystyle w_{ijj} = \exp \left( - \frac{1}{2} \Bigl( \frac{d_{ij}}{b} \Bigr)^2 \right)$

　 $d_{ij}$ は地点 $i, j$ の距離、 $b$ はバンド幅を表します。ただし、 $b$ は正規分布の標準偏差に対応します。それぞれ距離なので、0以上の値 $d_{ij} \geq 0$ 、また $b$ に関しては0除算防止のため、非負の値 $b \gt 0$ を満たす必要があります。
　ガウス型の重み関数は、0を除く0から1の値 $0 \lt w_{ijj} \leq 1$ をとります。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定
b = 40.0

# 対角要素を計算:(ガウス型)
w_diag = np.exp(-0.5 * (dist_vec / b)**2)
print(w_diag.round(2))
print(w_diag.shape)

# 重み行列を作成
W = np.diag(w_diag)
print(W.round(2))
print(W.shape)

[1.   0.27 0.79 0.1  0.65 1.   0.09 0.91 0.99 0.26]
(10,)
[[1.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.27 0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.79 0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.1  0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.65 0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   1.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.09 0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.91 0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.99 0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.26]]
(10, 10)

　バンド幅 $b$ を指定して、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。対角要素は $w_{ijj} = \exp(- \frac{1}{2} (\frac{d_{ij}}{b})^2)$ を計算します。非対角要素は $w_{ijl} = 0$ です。

bi-square型

　bi-square型の重み関数(bi-squareカーネル)による重み行列を作成します。

　bi-squareカーネルは、次の式で定義されます。

$\displaystyle w_{ijj} = \begin{cases} \Bigl( 1 - (\frac{d_{ij}}{b})^2 \Bigr)^2 &\quad (d_{ij} \lt b) \\ 0 &\quad (d_{ij} \geq b) \end{cases}$

　 $d_{ij}$ は地点 $i, j$ の距離、 $b$ はバンド幅を表します。それぞれ距離なので、0以上の値 $d_{ij} \geq 0$ 、また $b$ に関しては0除算防止のため、非負の値 $b \gt 0$ を満たす必要があります。
　bi-square型の重み関数は、0から1の値 $0 \leq w_{ijj} \leq 1$ をとります。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定
b = 40.0

# 対角要素を計算:(bi-square型)
w_diag = (1.0 - (dist_vec / b)**2)**2
print(w_diag.round(2))

# バンド幅外の要素を置換
w_diag[dist_vec >= b] = 0.0
print(w_diag.round(2))
print(w_diag.shape)

# 重み行列を作成
W = np.diag(w_diag)
print(W.round(2))
print(W.shape)

[ 0.98  2.51  0.27 13.31  0.02  0.99 14.62  0.66  0.95  2.93]
[0.98 0.   0.27 0.   0.02 0.99 0.   0.66 0.95 0.  ]
(10,)
[[0.98 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.27 0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.02 0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.99 0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.66 0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.95 0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]]
(10, 10)

　バンド幅 $b$ を指定して、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。対角要素は $w_{ijj} = (1 - (\frac{d_{ij}}{b})^2)^2$ を計算して、 $d_{ij} \geq b$ の要素を $w_{ijj} = 0$ に置き換えます。非対角要素は $w_{ijl} = 0$ です。

tri-cube型

　tri-cube型の重み関数(tri-cubeカーネル)を可視化します。

　tri-cubeカーネルは、次の式で定義されます。

$\displaystyle w_{ijj} = \begin{cases} \Bigl( 1 - (\frac{d_{ij}}{b})^3 \Bigr)^3 &\quad (d_{ij} \lt b) \\ 0 &\quad (d_{ij} \geq b) \end{cases}$

　 $d_{ij}$ は地点 $i, j$ の距離、 $b$ はバンド幅を表します。それぞれ距離なので、0以上の値 $d_{ij} \geq 0$ 、また $b$ に関しては0除算防止のため、非負の値 $b \gt 0$ を満たす必要があります。
　tri-cube型の重み関数は、0から1の値 $0 \leq w_{ijj} \leq 1$ をとります。

　バンド幅を指定して、重み行列を作成します。

# バンド幅を指定
b = 40.0

# 対角要素を計算:(tri-cube型)
w_diag = (1.0 - (dist_vec / b)**3)**3
print(w_diag.round(2))

# バンド幅外の要素を置換
w_diag[dist_vec >= b] = 0.0
print(w_diag.round(2))
print(w_diag.shape)

# 重み行列を作成
W = np.diag(w_diag)
print(W.round(2))
print(W.shape)

[ 1.0000e+00 -3.1500e+01  2.9000e-01 -7.3385e+02  1.0000e-02  1.0000e+00
 -8.8275e+02  7.7000e-01  9.9000e-01 -4.1570e+01]
[1.   0.   0.29 0.   0.01 1.   0.   0.77 0.99 0.  ]
(10,)
[[1.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.29 0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.01 0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   1.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.77 0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.99 0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]]
(10, 10)

　バンド幅 $b$ を指定して、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。対角要素は $w_{ijj} = (1 - (\frac{d_{ij}}{b})^3)^3$ を計算して、 $d_{ij} \geq b$ の要素を $w_{ijj} = 0$ に置き換えます。非対角要素は $w_{ijl} = 0$ です。

重み行列の実装

　次は、5つの重み関数の計算を行う関数を作成します。

　重み関数の計算と重み行列の作成を自作関数として定義します。

# 重み関数を定義
def get_weight_matrix(distance, bandwidth, function='bi-square', adjust_flag=False):

    # 変数を設定
    d = distance
    b = bandwidth

    # 対角要素を計算
    if function == 'moving window':

        # moving window型
        w = np.ones_like(d)
        w[d >= b] = 0.0

    elif function == 'exponential':

        # バンド幅を調整
        if adjust_flag:
            b *= 1.0 / 3.0

        # 指数型
        w = np.exp(- d / b)

    elif function == 'Gaussian':

        # バンド幅を調整
        if adjust_flag:
            b *= 1.0 / np.sqrt(6.0)

        # ガウス型
        w = np.exp(-0.5 * (d / b)**2)

    elif function == 'bi-square':

        # bi-square型
        w = (1.0 - (d / b)**2)**2
        w[d >= b] = 0.0

    elif function == 'tri-cube':

        # tri-cube型
        w = (1.0 - (d / b)**3)**3
        w[d >= b] = 0.0

    # 重み行列を作成
    W = np.diag(w)

    return W

　重み関数(カーネル)の計算処理を文字列 function で切り替えて、重み行列の対角要素 w を計算します。さらに、np.diag() で対角行列 W に変換して出力します。
　np.diag() に対角要素を渡すと対角行列を返し、対角行列を渡すと対角要素を返します。

　有効バンド幅を指定する(バンド幅の値を調整する)場合は、因子 adjust_flag で切り替えます。
　有効バンド幅については「GWRモデルの重み関数の可視化」を参照してください。

　実装した関数を使って、重み行列を作成します。

# 重み行列を作成
W = get_weight_matrix(distance=dist_vec, bandwidth=b, function='bi-square')
print(W.round(2))
print(W.shape)

[[0.98 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.27 0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.02 0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.99 0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.66 0.   0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.95 0.  ]
 [0.   0.   0.   0.   0.   0.   0.   0.   0.   0.  ]]
(10, 10)

　以上で、GWRモデルの重み行列の作成処理を実装できました。

　この記事では、GRMモデルの重み関数と重み行列をプログラムで確認しました。次からの記事では、図で確認します。

参考文献

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

作者:植杉　大
Independently published

Amazon

実践Data Scienceシリーズ Rではじめる地理空間データの統計解析入門

作者:村上大輔
講談社

Amazon

おわりに

　この先の記事で重み行列の計算が度々登場しそうだったので関数を切り替えられるように実装したのですが、本の内容的にはbi-square型しか使わないんですよね。というわけでこの先の記事ではこの自作関数は使いませんでした。アニメーションの作図には使ってます。
　この記事ごと要らないのではとも思ったのですが、可視化の方の記事で数式・プログラム・図の3つを扱うのは、書くのも読むのも重いよなと思い別記事として投稿しました。

　余談ですが、このブログの方向性が固まってからの記事では(概ね)、数式・アルゴリズム・プログラム・グラフのうち扱う内容によってタイトルが付けられています。「数式とプログラム」の対応を解説する記事では「○○の計算」さらにアルゴリズムも扱うと「○○の実装」、「数式とグラフ」の対応(式と図を見比べる目的)だと「○○の可視化」さらにプログラムも扱う(作図自体も目的)と「○○の作図」といった名付けルール(例外あり)でやってます。正直、ルールを優先して不自然な言い回しのタイトルになることもありますが、記事の管理を考えると構成などを共通化させたくなります。
　今回の記事だと、まぁ自作関数として実装しているからというのもありますが、「距離の計算→重みの計算→重み行列への変換」がアルゴリズムと(言うほどではないですが)みなしてる感じです。
　○○の憂鬱みたいなタイトルに心地よさを感じる世代なもんで。

【次の内容】

からっぽのしょこ

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

【Python】4.3.2-3：GWRモデルの重み行列の実装【空間データサイエンス入門のノート】

はじめに

4.3.2-3 GWRモデルの重み行列の実装

重み行列の数式表現

疑似データの生成

重み行列の計算

Moving Window型

指数型

ガウス型

bi-square型

tri-cube型

重み行列の実装

参考文献

おわりに