【Python】4.3.4-7：GWRモデルの最小二乗法の可視化【空間データサイエンス入門のノート】 (original) (raw)

はじめに

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

この記事では、GWRモデルの最小二乗法について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

• はじめに
• 4.3.4-7 GWRモデルの最小二乗法の可視化
  • 最小二乗法の図
    * 観測データの図
    * 重み行列の図
    * 係数パラメータの推定
• 参考文献
• おわりに

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

import numpy as np import libpysal as ps import geopandas as gpd import matplotlib.pyplot as plt

最小二乗法の図

まずは、GWRモデルに対するOLSの定義や性質を可視化します。

観測データの図

例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。

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

geom_df = gpd.read_file(ps.examples.get_path('G_utm.shp'))

geom_df = geom_df[['PctFB', 'PctBlack', 'PctRural', 'PctBach', 'X', 'Y', 'geometry']] geom_df

159 rows × 7 columns

libpysal ライブラリに含まれるジョージア州の空間データを利用します。
(利用する列を取り出す必要はないのですが、後ほど列を追加した際に確認しにくくなるので、不要な列を取り除いておきます。)

説明変数と被説明変数を設定します。

tmp_X = geom_df[['PctFB', 'PctBlack', 'PctRural']].values.astype(np.float64)

N, K = tmp_X.shape print(N, K)

X = np.hstack( [np.ones(shape=(N, 1)), tmp_X] ) print(X[:5].round(1)) print(X.shape)

y = geom_df['PctBach'].values.astype(np.float64)

K = X.shape[1] - 1

159 3
[[  1.    0.6  20.8  75.6]
 [  1.    1.6  26.9 100. ]
 [  1.    0.3  15.4  61.7]
 [  1.    0.1  51.7 100. ]
 [  1.    1.4  42.4  42.7]]
(159, 4)

Pandasデータフレームから各次元 の説明変数 として利用する列を取り出し、 番目の次元(定数項用)の値 と列方向に結合して、各地点 の説明変数 とします。また、 個の地点を行方向に結合して説明変数 とします。
同様に、 個の被説明変数 として用いる列を取り出して、NumPy配列に変換します。
データフレームの行数がデータ数 、取り出した列数がパラメータ数 に対応します。

説明変数と被説明変数のコロプレス図を作成します。

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

col_label_lt = ['PctBach', 'PctFB', 'PctBlack', 'PctRural']

R, C = 2, 2

fig, axes = plt.subplots(nrows=R, ncols=C, figsize=(11, 10), dpi=100, facecolor='White', constrained_layout=True) fig.suptitle('observation data', fontsize=20)

k = 0

for r in range(R): for c in range(C):

    col_label = col_label_lt[k]

    
    if k == 0:
        var_label = '$y_j$'
    else:
        var_label = f'$x_{{j{k}}}$'

    
    ax = axes[r, c]
    geom_df.plot(ax=ax, column=col_label, 
                 edgecolor='white', linewidth=0.5, 
                 legend=True, legend_kwds={'label': var_label}) 
    ax.set_xlabel('$u_j$')
    ax.set_ylabel('$v_j$')
    ax.set_title(col_label)
    ax.grid()
    ax.set_aspect('equal', adjustable='box')

    
    k += 1

plt.show()

変数のコロプレス図：元データ

個の地点について、被説明変数 と各次元の説明変数 の値をそれぞれコロプレス図で示します。

説明変数と被説明変数を標準化します。

tmp_X = (tmp_X - tmp_X.mean(axis=0)) / tmp_X.std(axis=0) y = (y - y.mean()) / y.std()

geom_df[['PctFB', 'PctBlack', 'PctRural']] = tmp_X geom_df['PctBach'] = y

X = np.hstack( [np.ones(shape=(N, 1)), tmp_X] ) print(X[:5].round(1)) print(X.mean(axis=0)) print(X.std(axis=0)) print(y[:5].round(1)) print(y.mean()) print(y.std())

[[ 1.  -0.4 -0.4  0.2]
 [ 1.   0.4 -0.   1.1]
 [ 1.  -0.7 -0.7 -0.3]
 [ 1.  -0.8  1.4  1.1]
 [ 1.   0.2  0.9 -1. ]]
[ 1.00000000e+00 -1.34064667e-16 -2.23441112e-17 -8.04388003e-16]
[0. 1. 1. 1.]
[-0.5 -0.8 -0.8 -0.3  0.4]
2.2344111187424534e-16
1.0

この例では作図(の配色の共通化)用に、 の次元(列)ごとと を標準化しておきます。ただし、0番目の次元(定数項用の変数)については、標準偏差が0のため、0除算になり計算できません。

標準化した説明変数と被説明変数をグラフで確認します。

変数のコロプレス図：標準化データ

の各次元(グラフ)と ごとに標本平均が0、標本標準偏差が1になるように標準化しました。各変数の値は変わりますが、大小関係などは変わりません。

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

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配列に格納します。
各地点 の経度を 、経度を として、位置を で表します。

ジョージア州における各地点のグラフを作成します。

fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', constrained_layout=True) geom_df.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5) ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], color='black', s=5) ax.set_xlabel('$u_j$') ax.set_ylabel('$v_j$') ax.set_title(f'$N = {N}$', loc='left') fig.suptitle('State of Georgia ', fontsize=20) ax.grid() ax.set_aspect('equal', adjustable='box') plt.show()

ジョージア州の各地点

区域ごとの点を各地域の位置として用います。

基準地点と各地点の距離の関係をアニメーションで確認します。

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

i = 0

fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', constrained_layout=True) fig.suptitle('State of Georgia ', fontsize=20)

def update(j):

ax.cla()


u_i, v_i = coord_arr[i]
u_j, v_j = coord_arr[j]


dist = np.sqrt(np.sum((coord_arr[j] - coord_arr[i])**2))


geom_df.plot(ax=ax, color='white', edgecolor='black', linewidth=0.5) 
ax.scatter(x=coord_arr[:, 0], y=coord_arr[:, 1], 
           color='black', s=5) 
ax.scatter(x=u_j, y=v_j, s=50, color='black', 
           label=f'$(u_j, v_j) = ({u_j:.1f}, {v_j:.1f})$') 
ax.scatter(x=u_i, y=v_i, s=50, color='red', 
           label=f'$(u_i, v_i) = ({u_i:.1f}, {v_i:.1f})$') 
ax.set_xlabel('$u_j$')
ax.set_ylabel('$v_j$')
ax.set_title(f'$N = {N}, i = {i+1}, j = {j+1}$', loc='left')
ax.grid()
ax.legend(loc='upper left')
ax.set_aspect('equal', adjustable='box')

ani = FuncAnimation(fig=fig, func=update, frames=N, interval=500)

ani.save( filename='georgia_dist.mp4', progress_callback = lambda i, n: print(f'frame: {i} / {n}') )

基準地点 を赤色の点、各地点 を黒色の点、2点間の距離 を青色の線分で示します。
GeoSeries.centroid.plot() は重心の点を描画するので、X, Y 列の点とはズレが生じます。

重み行列の図

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

b_idx = 116

b_lt = [] W_lt = []

for i in range(N):

dist_vec = np.sqrt(np.sum((coord_arr - coord_arr[i])**2, axis=1))


b = np.sort(dist_vec)[b_idx]


w_diag = (1.0 - (dist_vec / b)**2)**2
w_diag[dist_vec >= b] = 0.0


W_i = np.diag(w_diag)


b_lt.append(b)
W_lt.append(W_i)

print(W_lt[0].round(2)) print(W_lt[0].shape) print(len(W_lt))

[[1.   0.   0.   ... 0.   0.   0.  ]
 [0.   0.88 0.   ... 0.   0.   0.  ]
 [0.   0.   0.99 ... 0.   0.   0.  ]
 ...
 [0.   0.   0.   ... 0.07 0.   0.  ]
 [0.   0.   0.   ... 0.   0.53 0.  ]
 [0.   0.   0.   ... 0.   0.   0.52]]
(159, 159)
159

基準地点 ごとに、各地点 とのユークリッド距離 を計算し、重み行列 を作成して、リストに格納していきます。

全ての基準地点で共通のバンド幅 を指定します。この場合、基準地点ごとにバンド幅の範囲内の地点数が変わります。
または、バンド幅の範囲内に含まれる地点数 b_idx を指定して、基準地点ごとにバンド幅 を設定します。基準地点と各地点の距離を昇順に並べ替えて b_idx 番目の距離をバンド幅 として用います。これは mgwr ライブラリにおける処理と同じ操作です。
バンド幅を指定(固定)する場合は、コメントアウトしてください。

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

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

i = 0

W_i = W_lt[i]

geom_df['weight'] = np.diag(W_i)

b = b_lt[i]

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

fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', constrained_layout=True) geom_df.plot(ax=ax, column='weight', cmap='viridis', vmin=0.0, vmax=1.0, edgecolor='white', linewidth=0.5, legend=True, legend_kwds={'label': '$w_{ijj}$'}) 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('weight matrix: (bi-square kernel)', fontsize=20) fig.legend() 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()

重み行列(の対角要素)のコロプレス図

基準地点 の重み行列の対角要素 の値をコロプレス図で示します。
基準地点 の重みが で、離れた地点ほど値が0に近付きます。この例では、bi-square型の重み関数により重み付けしているので、バンド幅の範囲外は0になります。

重み付けした説明変数と被説明変数をグラフで確認します。

上段は元の変数 、下段は重み付けした変数 を表します。分かりやすいように、値が0のとき白色になるように調整しています。

基準地点ごとに各地点の重みが変わるので、各地点のモデルへの影響の度合いも変わります。

係数パラメータの推定

係数パラメータの推定値を計算します。

Beta_hat = np.zeros(shape=(N, K+1))

for i in range(N):

W_i = W_lt[i]


Beta_hat[i] = np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i @ y

print(Beta_hat[:5].round(2)) print(Beta_hat.shape)

[[-0.23  0.23  0.06 -0.43]
 [-0.28  0.17  0.1  -0.41]
 [-0.25  0.2   0.07 -0.43]
 [-0.23  0.15  0.05 -0.36]
 [ 0.19  0.72 -0.17 -0.24]]
(159, 4)

全ての基準地点をまとめた係数パラメータの推定値 として、基準地点 ごとに、係数パラメータの推定値 を で計算して、2次元配列に格納していきます。

係数パラメータの推定値のコロプレス図を作成します。

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

col_var_label_lt = ['intercept', 'PctFB', 'PctBlack', 'PctRural'] col_param_label_lt = ['param_itercept', 'param_PctFB', 'param_PctBlack', 'param_PctRural']

geom_df[col_param_label_lt] = Beta_hat

R, C = 2, 2

fig, axes = plt.subplots(nrows=R, ncols=C, figsize=(11, 10), dpi=100, facecolor='White', constrained_layout=True) fig.suptitle('estimated coefficient parameters', fontsize=20)

k = 0

for r in range(R): for c in range(C):

    col_var_label   = col_var_label_lt[k]
    col_param_label = col_param_label_lt[k]

    
    beta_max = max(Beta_hat[:, k].max(), abs(Beta_hat[:, k].min()))

    
    param_label = f'$\\beta_{{i{k}}}$'
    
    
    ax = axes[r, c]
    geom_df.plot(ax=ax, column=col_param_label, 
                 
                 edgecolor='white', linewidth=0.5, 
                 legend=True, legend_kwds={'label': param_label}) 
    ax.set_xlabel('$u_j$')
    ax.set_ylabel('$v_j$')
    ax.set_title(col_var_label)
    ax.grid()
    ax.set_aspect('equal', adjustable='box')

    
    k += 1

plt.show()

係数パラメータの推定値のコロプレス図

各地点(区域)が基準地点 に対応し、 個のグラフはそれぞれ各次元 の係数パラメータ の値をコロプレス図で示します。

被説明変数の推定値(理論値)を計算します。

Y_hat = X @ Beta_hat.T print(Y_hat.round(2)) print(Y_hat.shape)

[[-0.43 -0.46 -0.44 ... -0.08 -0.13 -0.44]
 [-0.62 -0.68 -0.65 ...  0.13  0.11 -0.62]
 [-0.3  -0.33 -0.31 ... -0.09 -0.16 -0.33]
 ...
 [-0.13 -0.11 -0.12 ... -0.24 -0.3  -0.13]
 [-0.81 -0.77 -0.8  ... -0.75 -0.72 -0.71]
 [-0.34 -0.35 -0.34 ... -0.15 -0.2  -0.34]]
(159, 159)

全ての基準地点をまとめた被説明変数の推定値 を で計算します。

被説明変数の推定値のコロプレス図を作成します。

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

geom_df['explained'] = np.diag(Y_hat)

fig, ax = plt.subplots(figsize=(10, 8), dpi=100, facecolor='white', constrained_layout=True) geom_df.plot(ax=ax, column='explained', edgecolor='white', linewidth=0.5, legend=True, legend_kwds={'label': '$\hat{y}_{ii}$'}) ax.set_xlabel('$u_j$') ax.set_ylabel('$v_j$') ax.set_title(f'$N = {N}$', loc='left') fig.suptitle('explained variables', fontsize=20) ax.grid() ax.set_aspect('equal', adjustable='box') plt.show()

基準地点の被説明変数の推定値のコロプレス図

各地点(区域)が基準地点 に対応し、 個の基準地点の被説明変数の推定値(理論値) の値をコロプレス図で示します。 の対角要素に対応します。

ここまでの記事では、GWRモデルの最小二乗法を確認しました。次の記事からはバンド幅の最適化を確認します。

参考文献

おわりに

間に合いませんでした。まだ作業中です。もう1・2枚図を追加するつもりでいます。

2024年10月1日は、えびちゅうこと私立恵比寿中学の桜井えまさんと仲村悠菜さんの加入2周年の日です。

お二人ともしっかり馴染んでるどころかキャラも役割もしっかり見えて頼もしい限りですね。

【次の内容】

つづく