からっぽのしょこ

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

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

はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

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()

重み行列のコロプレス図:moving window型

 基準地点  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()

重み行列のヒートマップ:moving window型

 同じ地点の組み合わせの要素(対角要素)  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 を赤色の破線で示します。

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

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

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

# ラベル用の文字列を指定
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()

重み行列のヒートマップ:ガウス型


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

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

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


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()

重み行列のコロプレス図:bi-square型

 基準地点  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()

重み行列のヒートマップ:bi-square型

 対角要素が  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()

重み行列のコロプレス図:tri-cube型

 基準地点  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})^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()

重み行列のヒートマップ:tri-cube型

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

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

 バンド幅  b が大きいほど、値を持つ重み  w_{ijj} \gt 0 の地点が増える(  w_{ijj} = 0 の地点が減る)ので、モデルに影響する地点が増え、また影響が強まります。

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


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

参考文献

おわりに

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

【次の内容】

www.anarchive-beta.com