はじめに

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

　この記事では、GWRモデルの重み行列について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
4.3.3 GWRモデルの重み行列の可視化
参考文献
おわりに

4.3.3 GWRモデルの重み行列の可視化

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

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

# ライブラリを読込
import numpy as np
import libpysal as ps
import geopandas as gpd
import matplotlib.pyplot as plt

重み行列の数式表現

　まずは、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$ です。

　この重み行列をコロプレス図やヒートマップで可視化します。

データの読込

　可視化の前に、例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。
　Georgiaデータセットについては「【Python】4.3.4-7：GWRモデルの最小二乗法の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

　Georgiaデータセットを読み込みます。

# 空間データを読込
geom_df = gpd.read_file(ps.examples.get_path('G_utm.shp'))
geom_df

	AREA	PERIMETER	G_UTM_	G_UTM_ID	Latitude	Longitud	TotPop90	PctRural	PctBach	PctEld	PctFB	PctPov	PctBlack	X	Y	AreaKey	geometry
0	1.331370e+09	207205.0	132	133	31.75339	-82.28558	15744	75.6	8.2	11.43	0.64	19.9	20.76	941396.6	3521764	13001	POLYGON ((931869.062 3545540.5, 934111.625 354...
1	8.929300e+08	154640.0	157	158	31.29486	-82.87474	6213	100.0	6.4	11.77	1.58	26.0	26.86	895553.0	3471916	13003	POLYGON ((867016.312 3482416, 884309.375 34826...
2	7.434020e+08	130431.0	148	146	31.55678	-82.45115	9566	61.7	6.6	11.11	0.27	24.1	15.42	930946.4	3502787	13005	POLYGON ((914656.875 3512190, 924718.375 35125...
3	9.053950e+08	185737.0	158	155	31.33084	-84.45401	3615	100.0	9.4	13.17	0.11	24.8	51.67	745398.6	3474765	13007	POLYGON ((744258.625 3480598.5, 765025.062 348...
4	6.941830e+08	151347.0	76	79	33.07193	-83.25085	39530	42.7	13.3	8.64	1.43	17.5	42.39	849431.3	3665553	13009	POLYGON ((832974.188 3677273.5, 834048.688 367...
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
154	7.528500e+08	174980.0	7	6	34.80497	-84.96616	72462	70.0	12.0	9.48	2.55	11.1	4.06	686891.4	3855274	13313	POLYGON ((684405.938 3873304.5, 684661.562 387...
155	9.935170e+08	160299.0	128	124	31.97034	-83.43574	7008	100.0	7.6	15.71	0.09	28.6	31.76	838551.5	3538547	13315	POLYGON ((845655.625 3557787.75, 846281.25 355...
156	1.230650e+09	180896.0	40	39	33.78664	-82.74436	10597	59.6	10.4	16.64	0.43	22.6	45.94	891228.5	3749769	13317	POLYGON ((902279.938 3768771.75, 906628.562 37...
157	1.175160e+09	161482.0	83	82	32.79853	-83.16759	10228	100.0	8.8	11.36	0.29	15.3	41.99	858796.9	3637891	13319	POLYGON ((867125.062 3651962.75, 867417.75 365...
158	1.489420e+09	200985.0	141	139	31.55269	-83.84816	19745	71.1	6.3	11.50	0.59	26.2	30.71	801018.1	3487328	13321	POLYGON ((798619.812 3525322, 799558.25 352473...

159 rows × 17 columns

　libpysal ライブラリに含まれるジョージア州の空間データを利用します。

　位置データを作成します。

# 位置情報を取得
coord_arr = np.stack(
    [geom_df['X'].values.astype(np.float64),
     geom_df['Y'].values.astype(np.float64)], 
    axis=1
)
print(coord_arr[:5])
print(coord_arr.shape)

[[ 941396.6 3521764. ]
 [ 895553.  3471916. ]
 [ 930946.4 3502787. ]
 [ 745398.6 3474765. ]
 [ 849431.3 3665553. ]]
(159, 2)

　Pandasデータフレームから経度 X 列と緯度 Y 列(に対応する値?これは何ですか?)を取り出してNumPy配列に格納します。
　各地点 $j$ の経度を $u_i$ 、経度を $v_i$ として、位置を $(u_i, v_i)$ で表します。

　基準地点を指定して、各地点までの距離を作成します。

# 地点数を取得
N = len(geom_df)

# 基準地点を指定
i = 0

# 距離を計算
dist_vec = np.sqrt(np.sum((coord_arr - coord_arr[i])**2, axis=1))
print(dist_vec[:5].round(1))
print(dist_vec.shape)

[     0.   67723.4  21664.1 201554.3 170683.6]
(159,)

　基準地点(の番号) $i$ を指定して、基準地点 $i$ と各地点 $j$ のユークリッド距離 $d_{ij} = \sqrt{(u_j - u_i)^2 + (v_j - v_i)^2}$ を計算します。

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

重み行列の図

　次は、GWRモデルで利用する重み関数の定義式を確認して、グラフを作成します。コロプレス図については本を参照してください。
　アニメーションの作図コードは「Introduction-to-SpatialDataScience/code/GWR/weight_matrix.py at main · anemptyarchive/Introduction-to-SpatialDataScience · GitHub」を参照してください。

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 = 270000.0

# 対角要素を計算:(moving window型)
w_diag = np.ones_like(dist_vec)
w_diag[dist_vec >= b] = 0.0

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

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

　バンド幅 $b$ を指定して、重み関数に従い対角要素 $w_{ijj}$ を計算し、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。
　np.diag() に対角要素を渡すと対角行列を返し、対角行列を渡すと対角要素を返します。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 重みを格納
geom_df['weight'] = np.diag(W)

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = 1 ~ (d_{ij} < b)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b)$'

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': fml_label}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, 
           color='red', s=50) # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('moving window weight function', fontsize=20)
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ を赤色の点、基準地点からのバンド幅 $b$ を赤色の破線で示します。

　バンド幅の範囲内 $d_{ij} \lt b$ の地点の重みは $w_{ijj} = 1$ 、範囲外 $d_{ij} \geq b$ の地点の重みは $w_{ijj} = 0$ になります。

　重み行列のヒートマップを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = 1 ~ (d_{ij} < b, j = l)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b, j = l)$\n'
fml_label += '$w_{ijl} = 0 ~ (j \\neq l)$'

# 重み行列を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
cs = ax.pcolor(W, vmin=0.0, vmax=1.0) # 行列
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$l$')
ax.set_ylabel('$j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (moving window)', fontsize=20)
fig.colorbar(cs, label=fml_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　同じ地点の組み合わせの要素(対角要素) $w_{ijj}$ のみが値をとるので、重み行列 $\mathbf{W}_i$ は対角行列になります。

　バンド幅と重み行列の関係をアニメーションで確認します。

　バンド幅 $b$ が大きいほど、重みが $w_{ijj} = 1$ の地点が増える( $w_{ijj} = 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 = 270000.0

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

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

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.78 0.   ... 0.   0.   0.  ]
 [0.   0.   0.92 ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.42 0.   0.  ]
 [0.   0.   0.   ... 0.   0.59 0.  ]
 [0.   0.   0.   ... 0.   0.   0.59]]
(159, 159)

　バンド幅 $b$ を指定して、重み関数に従い対角要素 $w_{ijj}$ を計算し、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 重みを格納
geom_df['weight'] = np.diag(W)

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# ラベル用の文字列を指定
fml_label = '$w_{ijj} = \\exp(- \\frac{d_{ij}}{b})$'

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': fml_label}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, color='red', s=50, 
           label='$(u_i, v_i)$') # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('exponential weight function', fontsize=20)
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ を赤色の点、基準地点からのバンド幅 $b$ を赤色の破線で示します。

　バンド幅の内外は直接影響せず、基準地点 $i$ に近い地点ほど重みが大きく( $w_{ijj} = 1$ に近く)、遠い地点ほど重みが小さく( $w_{ijj} = 0$ に近く)なります。
　バンド幅と有効範囲を対応させる方法は「GWRモデルの重み行列の実装」を参照してください。

　重み行列のヒートマップを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = \\exp(- \\frac{d_{ij}}{b}) ~ (j = l)$\n'
fml_label += '$w_{ijl} = 0 ~ (j \\neq l)$'

# 重み行列を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
cs = ax.pcolor(W, vmin=0.0, vmax=1.0) # 行列
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$l$')
ax.set_ylabel('$j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (exponential)', fontsize=20)
fig.colorbar(cs, label=fml_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　バンド幅と重み行列の関係をアニメーションで確認します。

　バンド幅 $b$ に応じて重み関数の曲線の形状が変わり、各地点の重み(重み行列の対角要素)の値が変化します。
　バンド幅 $b$ が大きいほど、各地点の重みが大きくなるので、モデルへの影響が強まります。

　地点と重み行列の関係をアニメーションで確認します。

ガウス型

　ガウス型の重み関数(ガウスカーネル・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 = 270000.0

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

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

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.97 0.   ... 0.   0.   0.  ]
 [0.   0.   1.   ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.69 0.   0.  ]
 [0.   0.   0.   ... 0.   0.87 0.  ]
 [0.   0.   0.   ... 0.   0.   0.87]]
(159, 159)

　バンド幅 $b$ を指定して、重み関数に従い対角要素 $w_{ijj}$ を計算し、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 重みを格納
geom_df['weight'] = np.diag(W)

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# ラベル用の文字列を指定
fml_label = '$w_{ijj} = \\exp(- \\frac{1}{2} (\\frac{d_{ij}}{b})^2)$'

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': fml_label}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, color='red', s=50, 
           label='$(u_i, v_i)$') # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('Gaussian weight function', fontsize=20)
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ を赤色の点、基準地点からのバンド幅 $b$ を赤色の破線で示します。

　重み行列のヒートマップを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = \\exp(- \\frac{1}{2} (\\frac{d_{ij}}{b})^2) ~ (j = l)$\n'
fml_label += '$w_{ijl} = 0 ~ (j \\neq l)$'

# 重み行列を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
cs = ax.pcolor(W, vmin=0.0, vmax=1.0) # 行列
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$l$')
ax.set_ylabel('$j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (Gaussian)', fontsize=20)
fig.colorbar(cs, label=fml_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　バンド幅と重み行列の関係をアニメーションで確認します。

　地点と重み行列の関係をアニメーションで確認します。

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 = 270000.0

# 対角要素を計算:(bi-square型)
w_diag = (1.0 - (dist_vec / b)**2)**2
w_diag[dist_vec >= b] = 0.0

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

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.88 0.   ... 0.   0.   0.  ]
 [0.   0.   0.99 ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.06 0.   0.  ]
 [0.   0.   0.   ... 0.   0.52 0.  ]
 [0.   0.   0.   ... 0.   0.   0.51]]
(159, 159)

　バンド幅 $b$ を指定して、重み関数に従い対角要素 $w_{ijj}$ を計算し、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 重みを格納
geom_df['weight'] = np.diag(W)

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = (1 - (\\frac{d_{ij}}{b})^2)^2 ~ (d_{ij} < b)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b)$'

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': fml_label}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, color='red', s=50, 
           label='$(u_i, v_i)$') # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('bi-square weight function', fontsize=20)
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ を赤色の点、基準地点からのバンド幅 $b$ を赤色の破線で示します。

　基準地点 $i$ に近い地点ほど重みが大きく( $w_{ijj} = 1$ に近く)、バンド幅 $b$ に近い地点ほど重みが小さく( $w_{ijj} = 0$ に近く)なり、バンド幅の範囲外 $d_{ij} \geq b$ の地点の重みは $w_{ijj} = 0$ です。

　重み行列のヒートマップを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = (1 - (\\frac{d_{ij}}{b})^2)^2 ~ (d_{ij} < b, j = l)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b, j = l)$\n'
fml_label += '$w_{ijl} = 0 ~ (j \\neq l)$'

# 重み行列を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
cs = ax.pcolor(W, vmin=0.0, vmax=1.0) # 行列
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$l$')
ax.set_ylabel('$j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (bi-square)', fontsize=20)
fig.colorbar(cs, label=fml_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　対角要素が $w_{ijj} = 0$ になる場合もあります。

　バンド幅と重み行列の関係をアニメーションで確認します。

　バンド幅 $b$ が大きいほど、値を持つ重み $w_{ijj} \gt 0$ の地点が増える( $w_{ijj} = 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 = 270000.0

# 対角要素を計算:(tri-cube型)
w_diag = (1.0 - (dist_vec / b)**3)**3
w_diag[dist_vec >= b] = 0.0

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

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.95 0.   ... 0.   0.   0.  ]
 [0.   0.   1.   ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.04 0.   0.  ]
 [0.   0.   0.   ... 0.   0.62 0.  ]
 [0.   0.   0.   ... 0.   0.   0.61]]
(159, 159)

　バンド幅 $b$ を指定して、重み関数に従い対角要素 $w_{ijj}$ を計算し、基準地点 $i$ の重み行列 $\mathbf{W}_i$ を作成します。

　重み行列の対角要素のコロプレス図を作成します。

作図コード(クリックで展開)

# 重みを格納
geom_df['weight'] = np.diag(W)

# 基準地点の位置を取得
u_i, v_i = coord_arr[i]

# バンド幅の範囲を計算
t_vec = np.linspace(start=0.0, stop=2.0*np.pi, num=361) # ラジアン
u_vec = u_i + b * np.cos(t_vec)
v_vec = v_i + b * np.sin(t_vec)

# グラフサイズを設定
u = 50000 # 切り下げ・切り上げの単位を指定
u_min = np.floor(geom_df.bounds.minx.min() /u)*u
u_max = np.ceil(geom_df.bounds.maxx.max()  /u)*u
v_min = np.floor(geom_df.bounds.miny.min() /u)*u
v_max = np.ceil(geom_df.bounds.maxy.max()  /u)*u

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = (1 - (\\frac{d_{ij}}{b})^3)^3 ~ (d_{ij} < b)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b)$'

# 重みのコロプレス図を作成
fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', 
                       constrained_layout=True)
geom_df.plot(ax=ax, column='weight', vmin=0.0, vmax=1.0, 
             edgecolor='white', linewidth=0.5, 
             legend=True, legend_kwds={'label': fml_label}) # 対角要素
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) # 各地点
ax.scatter(x=u_i, y=v_i, color='red', s=50, 
           label='$(u_i, v_i)$') # 基準地点
ax.plot(u_vec, v_vec, 
        color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('tri-cube weight function', fontsize=20)
ax.grid()
ax.set_xlim(xmin=u_min, xmax=u_max)
ax.set_ylim(ymin=v_min, ymax=v_max)
ax.set_aspect('equal', adjustable='box')
plt.show()

　基準地点 $i$ を赤色の点、基準地点からのバンド幅 $b$ を赤色の破線で示します。

　重み行列のヒートマップを作成します。

作図コード(クリックで展開)

# ラベル用の文字列を指定
fml_label  = '$w_{ijj} = (1 - (\\frac{d_{ij}}{b})^3)^3 ~ (d_{ij} < b, j = l)$\n'
fml_label += '$w_{ijj} = 0 ~ (d_{ij} \\geq b, j = l)$\n'
fml_label += '$w_{ijl} = 0 ~ (j \\neq l)$'

# 重み行列を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
cs = ax.pcolor(W, vmin=0.0, vmax=1.0) # 行列
ax.invert_yaxis() # 軸の反転
ax.set_xlabel('$l$')
ax.set_ylabel('$j$')
ax.set_title(f'$N = {N}, i = {i+1}, b = {b:.1f}$', loc='left')
fig.suptitle('weight matrix: (tri-cube)', fontsize=20)
fig.colorbar(cs, label=fml_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　対角要素が $w_{ijj} = 0$ になる場合もあります。

　バンド幅と重み行列の関係をアニメーションで確認します。

　地点と重み行列の関係をアニメーションで確認します。

　ここまでの記事では、GWRモデルの重み行列を確認しました。次からの記事では、GWRモデルを確認していきます。

参考文献

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

作者:植杉　大
Independently published

Amazon

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

作者:村上大輔
講談社

Amazon

おわりに

　GWRで使う重み行列は、2章で扱った隣接関係による重み行列よりもシンプルなんですが、モデルに対してどう影響しているのかをイマイチ掴めなかったので、とりあえず図を描いてみました。
　行列と言っても意味を持つのは対角要素のみなので、扱い方は1次元的です。行列にしておくとモデルを行列計算で表現(処理)できてなんだろうなと思います。
　重い内容ではないので重み行列に関してはこの記事だけのつもりだったのですが、関数が5つに増えたあたりで重み関数の内容を別記事に切り離し、パラメータ推定の記事との兼ね合いで実装の内容を切り離しで、3記事になりました。

【次の内容】

からっぽのしょこ

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

【Python】4.3.3：GWRモデルの重み行列の可視化【空間データサイエンス入門のノート】

はじめに

4.3.3 GWRモデルの重み行列の可視化

重み行列の数式表現

データの読込

重み行列の図

Moving Window型

指数型

ガウス型

bi-square型

tri-cube型

参考文献

おわりに