はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
4.3.2 GWRモデルの重み関数の可視化
- 重み関数の図
- 重み関数の比較
参考文献
おわりに

4.3.2 GWRモデルの重み関数の可視化

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

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

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

重み関数の図

　まずは、GWRモデルで利用する重み関数の定義式を確認して、グラフを作成します。

Moving Window型

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

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

$\displaystyle f(d_{ij}) = \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値 $f(d_{ij}) \in {0, 1}$ をとります。

　moving windowカーネルをグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(moving window型)
weight_line = np.ones_like(dist_line)
weight_line[dist_line >= b] = 0.0

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

# ラベル用の文字列を指定
fml_label  = '$w = 1 ~ (d < b)$\n'
fml_label += '$w = 0 ~ (d \\geq b)$'

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('moving window weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を赤色の破線で示します。

　基準地点 $i$ と各地点 $j$ の距離 $d_{ij}$ がバンド幅 $b$ 未満 $d_{ij} \lt b$ であれば重みは $w_{ij} = 1$ 、バンド幅以上 $d_{ij} \geq b$ であれば $w_{ij} = 0$ になります。

指数型

　指数型の重み関数(指数カーネル・exponential kernel)を可視化します。

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

$\displaystyle f(d_{ij}) = \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 f(d_{ij}) \leq 1$ をとります。

　指数カーネルをグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(指数型)
weight_line = np.exp(- dist_line / b)

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

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

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('exponential weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を赤色の破線で示します。
　基準地点 $i$ と各地点 $j$ の距離 $d_{ij}$ が大きくなる(離れる)ほど重み $w_{ij}$ が小さくなります。基準の地点 $i$ (距離が $d_{ii} = 0$ )のとき $w_{ii} = \exp(0) = 1$ 、地点 $j$ の距離がバンド幅 $d_{ij} = b$ のとき $w_{ij} = \exp(- 1) = 0.367\dots$ となります。
　距離 $d_{ij}$ がどれだけ大きくなっても重みは $0$ にはなりません。

　有効バンド幅を指定した指数カーネルをグラフで確認します。

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

# 有効バンド幅を指定
b_star = 40.0

# バンド幅を計算
b = 1.0/3.0 * b_star

# 重み関数を計算:(指数型)
weight_line = np.exp(- dist_line / b)

# ラベル用の文字列を指定
param_label  = '$b = \\frac{1}{3} b^{*} = ' + f'{b:.1f}, '
param_label += 'b^{*} = ' + f'{b_star:.1f}$'

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=[b, b_star], ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle=['dashed', 'dotted']) # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(param_label, loc='left')
fig.suptitle('exponential weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を破線、有効バンド $b^{*} = 3 b$ を点線で示します。

　有効バンド幅 $b^{*}$ を指定し、バンド幅 $b = \frac{1}{3} b^{*}$ に変換して、重みを計算することでモデルに影響する範囲を調整できます。

ガウス型

　ガウス型の重み関数(ガウスカーネル・Gaussian kernel)を可視化します。

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

$\displaystyle f(d_{ij}) = \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 f(d_{ij}) \leq 1$ をとります。

　ガウスカーネルをグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(ガウス型)
weight_line = np.exp(-0.5 * (dist_line / b)**2)

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

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

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('Gaussian weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を赤色の破線で示します。
　距離は0以上の値 $d_{ij} \geq 0$ なので、平均 $\mu = 0$ 、標準偏差 $\sigma = b$ の正規化項の無いガウス分布(正規分布)の正の範囲と言えます。

　基準地点 $i$ と各地点 $j$ の距離 $d_{ij}$ が大きくなる(離れる)ほど重み $w_{ij}$ が小さくなります。基準の地点 $i$ (距離が $d_{ii} = 0$ )のとき $w_{ii} = \exp(0) = 1$ 、地点 $j$ の距離がバンド幅 $d_{ij} = b$ のとき $w_{ij} = \exp(- \frac{1}{2}) = 0.606\dots$ となります。
　距離 $d_{ij}$ がどれだけ大きくなっても重みは $0$ にはなりません。

　有効バンド幅を指定したガウスカーネルをグラフで確認します。

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

# 有効バンド幅を指定
b_star = 40.0

# バンド幅を計算
b = 1.0/np.sqrt(6.0) * b_star

# 重み関数を計算:(ガウス型)
weight_line = np.exp(-0.5 * (dist_line / b)**2)

# ラベル用の文字列を指定
param_label  = '$b = \\frac{1}{\\sqrt{6}} b^{*} = ' + f'{b:.1f}, '
param_label += 'b^{*} = ' + f'{b_star:.1f}$'

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=[b, b_star], ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle=['dashed', 'dotted']) # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(param_label, loc='left')
fig.suptitle('Gaussian weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を破線、有効バンド $b^{*} = \sqrt{6} b$ を点線で示します。

　有効バンド幅 $b^{*}$ を指定し、バンド幅 $b = \frac{1}{\sqrt{6}} b^{*}$ に変換して、重みを計算することでモデルに影響する範囲を調整できます。

bi-square型

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

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

$\displaystyle f(d_{ij}) = \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 f(d_{ij}) \leq 1$ をとります。

　bi-squareカーネルをグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(bi-square型)
tmp_weight_line = (1.0 - (dist_line / b)**2)**2
weight_line = np.zeros_like(dist_line)
weight_line[dist_line < b] = tmp_weight_line[dist_line < b]

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

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

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線:(バンド幅あり)
ax.plot(dist_line, tmp_weight_line, 
        color='black', linestyle='dotted', label=fnc_label) # 重み曲線:(バンド幅なし)
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('bi-square weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を赤色の破線で、また目安として関数 $f(d) = (1 - (\frac{d}{b})^2)^2$ の曲線を黒色の点線で示します。

　基準の地点 $i$ (距離が $d_{ii} = 0$ )のとき $w_{ii} = (1 - 0)^2 = 1$ 、地点 $j$ の距離がバンド幅 $d_{ij} = b$ のとき $w_{ij} = (1 - 1)^2 = 0$ となります。また、バンド幅以上 $d_{ij} \geq b$ のとき $w_{ij} = 0$ です。

tri-cube型

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

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

$\displaystyle f(d_{ij}) = \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 f(d_{ij}) \leq 1$ をとります。

　tri-cubeカーネルをグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(tri-cube型)
tmp_weight_line = (1.0 - (dist_line / b)**3)**3
weight_line = np.zeros_like(dist_line)
weight_line[dist_line < b] = tmp_weight_line[dist_line < b]

# グラフサイズを設定
w_min = -1.1
w_max = 1.1

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

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 10), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
ax.plot(dist_line, weight_line, 
        color='C0', label=fml_label) # 重み曲線:(バンド幅あり)
ax.plot(dist_line, tmp_weight_line, 
        color='black', linestyle='dotted', label=fnc_label) # 重み曲線:(バンド幅なし)
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('tri-cube weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　バンド幅 $b$ を赤色の破線で、また目安として関数 $f(d) = (1 - (\frac{d}{b})^3)^3$ の曲線を黒色の点線で示します。

　基準の地点 $i$ (距離が $d_{ii} = 0$ )のとき $w_{ii} = (1 - 0)^3 = 1$ 、地点 $j$ の距離がバンド幅 $d_{ij} = b$ のとき $w_{ij} = (1 - 1)^3 = 0$ となります。また、バンド幅以上 $d_{ij} \geq b$ のとき $w_{ij} = 0$ です。

重み関数の比較

　次は、5つの重み関数の曲線を並べて比較します。
　アニメーションの作図コードは「Introduction-to-SpatialDataScience/code/GWR/weight_function.py at main · anemptyarchive/Introduction-to-SpatialDataScience · GitHub」を参照してください。

　重み関数をグラフで確認します。

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

# バンド幅を指定
b = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算
weight_line_dic = {
    'moving window': np.ones_like(dist_line), 
    'exponential':   np.exp(- dist_line / b), 
    'Gaussian':      np.exp(-0.5 * (dist_line / b)**2), 
    'bi-square':     (1.0 - (dist_line / b)**2)**2, 
    'tri-cube':      (1.0 - (dist_line / b)**3)**3
}
weight_line_dic['moving window'][dist_line >= b] = 0.0
weight_line_dic['bi-square'][dist_line >= b] = 0.0
weight_line_dic['tri-cube'][dist_line >= b] = 0.0

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b_star, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
for fnc_label, weight_line in weight_line_dic.items():
    ax.plot(dist_line, weight_line, 
            label=fnc_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b = {b:.1f}$', loc='left')
fig.suptitle('weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　関数ごとに減衰の様子が異なります。

　バンド幅と曲線の形状の関係をアニメーションで確認します。

　バンド幅が大きいほど、基準の地点から遠い地点の重みが大きくなり、広い範囲の影響を受けます。

　有効バンド幅を指定した重み関数をグラフで確認します。

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

# 有効バンド幅を指定
b_star = 40.0

# 距離の範囲を指定
dist_line = np.linspace(start=0.0, stop=100.0, num=1001)

# 重み関数を計算:(有効バンド幅を調整)
weight_line_dic = {
    'moving window': np.ones_like(dist_line), 
    'exponential':   np.exp(- dist_line * 3.0/b_star), 
    'Gaussian':      np.exp(-3.0 * (dist_line / b_star)**2), 
    'bi-square':     (1.0 - (dist_line / b_star)**2)**2, 
    'tri-cube':      (1.0 - (dist_line / b_star)**3)**3
}
weight_line_dic['moving window'][dist_line >= b_star] = 0.0
weight_line_dic['bi-square'][dist_line >= b_star] = 0.0
weight_line_dic['tri-cube'][dist_line >= b_star] = 0.0

# グラフサイズを設定
w_min = -0.1
w_max = 1.1

# 重み関数を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', 
                       constrained_layout=True)
ax.vlines(x=b, ymin=w_min, ymax=w_max, 
          color='red', linewidth=2.0, linestyle='dashed') # バンド幅
for fnc_label, weight_line in weight_line_dic.items():
    ax.plot(dist_line, weight_line, 
            label=fnc_label) # 重み曲線
ax.set_xlabel('distance ($d$)')
ax.set_ylabel('weight ($w$)')
ax.set_title(f'$b^{{*}} = {b_star:.1f}$', loc='left')
fig.suptitle('weight function', fontsize=20)
ax.set_ylim(ymin=w_min, ymax=w_max)
ax.grid()
ax.legend()
plt.show()

　有効バンド幅と曲線の形状の関係をアニメーションで確認します。

　この記事では、GWRモデルの重み関数を図で確認しました。次の記事では、重み行列を図で確認します。

参考文献

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

作者:植杉　大
Independently published

Amazon

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

作者:村上大輔
講談社

Amazon

おわりに

　いくつか種類があるよと言われると全部試してみたくなりますよね。式だけ見せられると図にしてみたくなりますよね。というわけで重み関数をグラフにしてみました。ちなみに本に載ってるのは3つで、使うのは1つです(4章現在)。

　GWRのパラメータ推定の式変形でどうにも分からない部分があったので別の本を少し眺めてみたところ、他に2つ登場したのでそれもやりました。
　正直、重み行列については書き終わった後に読んだので見なかったことにしようかとも思ったのですが、既に描いてた図の軸範囲やバンド幅が全く同じだったことを縁に感じて追加して、重み関数のみを扱う記事に独立しました。
　まぁ良い図を求めて組み合わせを試行錯誤したところ、範囲100に対してバンド幅40が曲線の見栄え的にも性質の説明的にもいい感じなのは分かります。

　目的の情報は得られませんでしたが、、GWRに関しては数式はなく推定や可視化のRコードが充実していました(他の内容については読んでません)。

　2024年9月25日は、ukkaの結城りなさんと葵るりさんの加入3周年の日です。

　さらにニューシングルのリリース日でもあるので、その中から重い人の曲を最後に貼っておきます。
　私は、(えびちゅうを経由して)結城りなさんきっかけでukkaを追い始めて1年目のド新参です。よろしくお願いします！

【次の内容】

からっぽのしょこ

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

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

はじめに

4.3.2 GWRモデルの重み関数の可視化

重み関数の図

Moving Window型

指数型

ガウス型

bi-square型

tri-cube型

重み関数の比較

参考文献

おわりに