はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
4.3.4-7 GWRモデルの最小二乗法の実装
参考文献
おわりに

4.3.4-7 GWRモデルの最小二乗法の実装

　地理的加重回帰モデル(GWRモデル・geographically weighted regression model)に対する最小二乗法(OLS・ordinary least square)におけるパラメータの推定処理を実装します。
　地理的加重回帰モデルの定義や計算式については「4.3.4-7：GWRモデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」、線形回帰モデルや最小二乗法については「【Python】3.2.1-3：線形回帰モデルの最小二乗法の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

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

データの読込

　可視化の前に、例として利用する空間データを読み込んで前処理を行います。詳しくは本を参照してください。
　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 ライブラリに含まれるジョージア州の空間データを利用します。

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

# データ数を取得
N = len(geom_df)
print(N)

# 説明変数を取得
X = np.hstack(
    [np.ones(shape=(N, 1)), 
     geom_df[['PctFB', 'PctBlack', 'PctRural']].values.astype(np.float64)]
)
print(X[:5].round(1))
print(X.shape)

# 被説明変数を取得
y = geom_df['PctBach'].values.astype(np.float64).reshape((N, 1))
print(y[:5].round(1))
print(y.shape)

# パラメータ数を取得
K = X.shape[1] - 1
print(K)

159
[[  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)
[[ 8.2]
 [ 6.4]
 [ 6.6]
 [ 9.4]
 [13.3]]
(159, 1)
3

　Pandasデータフレームから各次元 $k$ の説明変数 $x_{1k}, \cdots, x_{Nk}$ として利用する列を取り出し、 $0$ 番目の次元(定数項用)の値 $x_{j0} = 1$ と列方向に結合して、各地点 $j$ の説明変数 $\mathbf{x}_j = (x_{j0}, x_{j1}, \cdots, x_{jK})$ とします。また、 $N$ 個の地点を行方向に結合して説明変数 $\mathbf{X} = (x_{10}, \cdots, x_{NK})$ とします。
　同様に、 $N$ 個の被説明変数 $\mathbf{y} = (y_1, \cdots, y_N)$ として用いる列を取り出して、NumPy配列に変換します。
　データフレームの行数がデータ数 $N$ 、取り出した列数がパラメータ数 $K$ に対応します。

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

# 位置情報を取得
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)$ で表します。

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

最小二乗法の計算

　次は、GWRモデルに対するOLSによるパラメータ推定の処理を確認します。

重み行列の設定

　重み行列については「【Python】4.3.3：GWRモデルの重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」、重み行列の作成については「【Python】4.3.2-3：GWRモデルの重み行列の実装【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

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

# バンド幅を指定:(バンド幅を固定する場合)
#b = 270000.0

# バンド幅内の地点数を指定:(地点数を固定する場合)
b_idx = 116

# 重み行列の受け皿を初期化
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]
    
    # 対角要素を計算:(bi-square型)
    w_diag = (1.0 - (dist_vec / b)**2)**2
    w_diag[dist_vec >= b] = 0.0
    
    # 重み行列を作成
    W_i = np.diag(w_diag)

    # 重み行列を格納
    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

　基準地点 $i$ ごとに、各地点 $j$ とのユークリッド距離 $d_{ij} = \sqrt{(u_j - u_i)^2 + (v_j - v_i)^2}$ を計算し、重み行列 $\mathbf{W}_i = (w_{i11}, \cdots, w_{iNN})$ を作成して、リストに格納していきます。

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

係数パラメータの推定

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

# 係数パラメータの受け皿を初期化
beta_hat_lt = []

# 基準地点ごとに処理
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

    # 係数パラメータを格納
    beta_hat_lt.append(beta_hat_i)
print(beta_hat_lt[0].round(2))
print(beta_hat_lt[0].shape)
print(len(beta_hat_lt))

[[14.22]
 [ 1.05]
 [ 0.02]
 [-0.09]]
(4, 1)
159

　基準地点 $i$ ごとに、係数パラメータの推定値 $\hat{\boldsymbol{\beta}}_i = (\hat{\beta}_{i0}, \hat{\beta}_{i1}, \cdots, \hat{\beta}_{iK})$ を次の式で計算して、リストに格納していきます。

$\displaystyle \hat{\boldsymbol{\beta}}_i = (\mathbf{X}^{\top} \mathbf{W}_i \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{W}_i \mathbf{y}$

　逆行列は linalg モジュールの inv() で計算できます。

　次の2つの処理は、計算処理の確認のために行います。計算結果は、パラメータ推定には使いません。

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

# 被説明変数の推定値の受け皿を初期化
y_hat_lt = []

# 基準地点ごとに処理
for i in range(N):

    # 係数パラメータの推定値を取得
    beta_hat_i = beta_hat_lt[i]
    
    # 被説明変数の推定値を計算
    y_hat_i = X @ beta_hat_i

    # 被説明変数の推定値を格納
    y_hat_lt.append(y_hat_i)
print(y_hat_lt[0][:5].round(2))
print(y_hat_lt[0].shape)
print(len(y_hat_lt))

[[ 8.5 ]
 [ 7.42]
 [ 9.26]
 [ 6.34]
 [12.69]]
(159, 1)
159

　基準地点 $i$ ごとに、被説明変数の推定値 $\hat{\mathbf{y}}_i = (\hat{y}_{i1}, \cdots, \hat{y}_{iN})$ を次の式で計算して、リストに格納していきます。

$\displaystyle \hat{\mathbf{y}}_i = \mathbf{X} \hat{\boldsymbol{\beta}}_i$

　基準地点ごとによる重み付け残差を計算します。

# 残差の受け皿を初期化
e_i_lt = []

# 基準地点ごとに処理
for i in range(N):

    # 重み行列を取得
    W_i = W_lt[i]
    
    # 被説明変数の推定値を取得
    y_hat_i = y_hat_lt[i]
    
    # 残差を計算
    e_i = W_i @ (y - y_hat_i)

    # 残差を格納
    e_i_lt.append(e_i)
print(e_i_lt[0][:5].round(2))
print(e_i_lt[0].shape)
print(len(e_i_lt))

[[-0.3 ]
 [-0.9 ]
 [-2.63]
 [ 0.63]
 [ 0.23]]
(159, 1)
159

　基準地点 $i$ ごとに、残差 $\mathbf{e}_i = (e_{i1}, \cdots, e_{iN})$ を次の式で計算して、リストに格納していきます。

$\displaystyle \mathbf{e}_i = \mathbf{W}_i (\mathbf{y} - \hat{\mathbf{y}}_i)$

分散パラメータの推定

　ハット行列を計算します。

# ハット行列を初期化
R = np.zeros(shape=(N, N))

# 基準地点ごとに処理
for i in range(N):

    # 重み行列を取得
    W_i = W_lt[i]

    # ハット行列の行ベクトルを計算
    r_i = X[i, np.newaxis] @ np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i
    
    # 行ベクトルを格納
    R[i] = r_i.flatten()
print(R.round(2))
print(R.shape)

[[ 0.04  0.02  0.05 ...  0.    0.01  0.01]
 [ 0.02  0.08 -0.   ... -0.    0.02  0.01]
 [ 0.05 -0.01  0.08 ...  0.   -0.    0.02]
 ...
 [ 0.    0.    0.   ...  0.07  0.01  0.  ]
 [ 0.    0.    0.   ...  0.01  0.06  0.  ]
 [ 0.02  0.01  0.02 ...  0.    0.01  0.03]]
(159, 159)

　ハット行列 $\mathbf{R} = (r_{11}, \cdots, r_{NN})$ を作成します。
　基準地点 $i$ ごとに、各行 $\mathbf{r}_i = (r_{i1}, \cdots, r_{iN})$ を次の式で計算して、配列に格納していきます。

$\displaystyle \begin{aligned} \mathbf{r}_i^{\top} &= \mathbf{x}_i^{\top} (\mathbf{X}^{\top} \mathbf{W}_i \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{W}_i \\ \mathbf{R} &= \begin{pmatrix} \mathbf{r}_1^{\top} \\ \mathbf{r}_2^{\top} \\ \vdots \\ \mathbf{r}_N^{\top} \end{pmatrix} \end{aligned}$

　基準地点ごとの(重み付けしていない)残差を計算します。

# 残差を計算
e = (np.identity(N) - R) @ y
print(e[:5].round(2))
print(e.shape)

[[-0.3 ]
 [-0.7 ]
 [-2.6 ]
 [ 2.33]
 [-0.28]]
(159, 1)

　基準地点 $i$ ごとの残差( $N$ 次元ベクトル) $\mathbf{e} = (e_{11}, \cdots, e_{NN})$ を次の式で計算します。

$\displaystyle \mathbf{e} = (\mathbf{I} - \mathbf{R}) \mathbf{y}$

　単位行列は np.identity() で作成できます。

　有効パラメータ数を計算します。

# 有効パラメータ数を計算
v1 = np.trace(R)
v2 = np.trace(R.T @ R)
print(v1.round(2))
print(v2.round(2))

11.8
8.29

　有効パラメータ数 $2 v_1 - v_2$ に関する値 $v_1, v_2$ を次の式で計算します。

$\displaystyle \begin{aligned} v_1 &= \mathrm{trace}(\mathbf{R}) \\ v_2 &= \mathrm{trace}(\mathbf{R}^{\top} \mathbf{R}) \end{aligned}$

　分散パラメータの推定値を作成します。

# 残差平方和を計算
rss = e.T @ e
print(rss.round(2))

# 分散パラメータの推定値を計算
sigma2_hat = rss / (N - (2.0*v1 - v2))
print(sigma2_hat.round(2))

# 標準偏差パラメータの推定値を計算
sigma_hat = np.sqrt(sigma2_hat)
print(sigma_hat.round(2))

[[1650.86]]
[[11.49]]
[[3.39]]

　分散パラメータの推定値 $\hat{\sigma}^2$ と標準偏差パラメータの推定値 $\hat{\sigma}$ を次の式で計算します。

$\displaystyle \begin{aligned} \hat{\sigma}^2 &= \frac{ \mathbf{e}^{\top} \mathbf{e} }{ N - (2 v_1 - v_2) } \\ \hat{\sigma} &= \sqrt{\hat{\sigma}^2} \end{aligned}$

係数パラメータの分散共分散行列の推定

　係数パラメータの推定値の分散共分散行列を作成します。

# 分散共分散行列の受け皿を初期化
Sigma_beta_hat_lt = []

# 基準地点ごとに処理
for i in range(N):
    
    # 重み行列を取得
    W_i = W_lt[i]
    
    # 中間変数を計算
    C_i = np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i
    
    # 係数パラメータの推定値の分散共分散行列を計算
    Sigma_beta_hat_i = sigma2_hat * C_i @ C_i.T

    # 分散共分散行列を格納
    Sigma_beta_hat_lt.append(Sigma_beta_hat_i)
print(Sigma_beta_hat_lt[0].round(3))
print(Sigma_beta_hat_lt[0].shape)
print(len(Sigma_beta_hat_lt))

[[ 3.61e+00 -6.31e-01 -3.40e-02 -2.50e-02]
 [-6.31e-01  2.71e-01  3.00e-03  4.00e-03]
 [-3.40e-02  3.00e-03  1.00e-03  0.00e+00]
 [-2.50e-02  4.00e-03  0.00e+00  0.00e+00]]
(4, 4)
159

　基準地点 $i$ ごとに、係数パラメータの推定値 $\hat{\boldsymbol{\beta}}_i$ の分散共分散行列を次の式で計算して、リストに格納していきます。

$\displaystyle \begin{aligned} \mathbf{C}_i &= (\mathbf{X}^{\top} \mathbf{W}_i \mathbf{X})^{-1} \mathbf{X}^{\top} \mathbf{W}_i \\ \mathbb{E} \Bigl[ (\hat{\boldsymbol{\beta}}_i - \mathbb{E}[\hat{\boldsymbol{\beta}}_i]) (\hat{\boldsymbol{\beta}}_i - \mathbb{E}[\hat{\boldsymbol{\beta}}_i])^{\top} \Bigr] &= \sigma^2 \mathbf{C}_i \mathbf{C}_i^{\top} \end{aligned}$

　以上の処理を実装します。

最小二乗法の実装

　続いて、GWRモデルに対するOLSによるパラメータ推定を行うクラスを作成します。

　OLSによるパラメータ推定を自作クラスとして定義します。

# GWRモデルに対するOLSを実装
class GWR:

    # データの設定
    def __init__(self, X, y, W):
        
        # 観測データを格納
        self.X = X # 説明変数
        self.y = y # 被説明変数
        self.W = W # 重み行列
        self.N = X.shape[0] # 観測データ数
        self.K = X.shape[1] - 1 # 係数パラメータ数

    # 係数パラメータの推定
    def get_beta_hat(self):
        
        # 観測データを取得
        X    = self.X
        y    = self.y.flatten() # (横ベクトルで扱う場合)
        W_lt = self.W
        
        # 基準地点ごとに処理
        beta_hat_lt = []
        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
        
            # 係数パラメータを格納
            beta_hat_lt.append(beta_hat_i)
        
        return beta_hat_lt
    
    # 誤差項の分散パラメータの推定
    def get_sigma2_hat(self):
        
        # 観測データを取得
        X    = self.X
        y    = self.y.flatten() # (横ベクトルで扱う場合)
        W_lt = self.W
        
        # 基準地点ごとに処理
        R = np.zeros(shape=(N, N)) # ハット行列を初期化
        for i in range(N):
        
            # 重み行列を取得
            W_i = W_lt[i]
        
            # ハット行列の行ベクトルを計算
            r_i = X[i, np.newaxis] @ np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i
            
            # 行ベクトルを格納
            R[i] = r_i.flatten()

        # 有効パラメータ数を計算
        v1 = np.trace(R)
        v2 = np.trace(R.T @ R)
        
        # 残差を計算
        e = (np.identity(N) - R) @ y
        
        # 残差平方和を計算
        rss = e.T @ e
        
        # 分散パラメータの推定値を計算
        sigma2_hat = rss / (self.N - (2.0*v1 - v2))
        
        return sigma2_hat
    
    # 係数パラメータの推定値の分散共分散行列の推定
    def get_sigma_mat_hat(self):

        # 観測データを取得
        X    = self.X
        W_lt = self.W

        # 分散パラメータの推定値を計算
        sigma2_hat = self.get_sigma2_hat()
        
        # 基準地点ごとに処理
        Sigma_beta_hat_lt = []
        for i in range(N):
            
            # 重み行列を取得
            W_i = W_lt[i]
            
            # 中間変数を計算
            C_i = np.linalg.inv(X.T @ W_i @ X) @ X.T @ W_i
            
            # 係数パラメータの推定値の分散共分散行列を計算
            Sigma_beta_hat_i = sigma2_hat * C_i @ C_i.T
        
            # 分散共分散行列を格納
            Sigma_beta_hat_lt.append(Sigma_beta_hat_i)

        return Sigma_beta_hat_lt

　出力変数 y が (N, 1) の2次元配列(縦ベクトル)の場合、.flatten() で (N,) の1次元配列(横ベクトル)に変形して処理します。y が (N,) の1次元配列の場合はそのまま処理されます。
　y を縦ベクトルで扱う( .flatten() を使わない)場合は、(残差平方和 rss や)分散パラメータの推定値 sigma2_hat が (1, 1) の2次元配列になるので、.item() でスカラに変換して出力する処理を追加する必要があります。

　重み行列はクラスの外で作成してインスタンス作成時に引数に渡すように実装しています。

　実装したクラスを使って、OLSによるパラメータ推定の計算をします。

# GWRのインスタンスを作成
gwr = GWR(X, y, W_lt)

# 係数パラメータの推定値を計算
beta_hat_lt = gwr.get_beta_hat()
print(beta_hat_lt[0].round(2))
print(beta_hat_lt[0].shape)
print(len(beta_hat_lt))

# 分散パラメータの推定値を計算
sigma2_hat = gwr.get_sigma2_hat()
print(sigma2_hat.round(2))

# 係数パラメータの推定値の分散共分散行列を計算
Sigma_beta_hat_lt = gwr.get_sigma_mat_hat()
print(Sigma_beta_hat_lt[0].round(3))
print(Sigma_beta_hat_lt[0].shape)
print(len(Sigma_beta_hat_lt))

[14.22  1.05  0.02 -0.09]
(4,)
159
11.49
[[ 3.61e+00 -6.31e-01 -3.40e-02 -2.50e-02]
 [-6.31e-01  2.71e-01  3.00e-03  4.00e-03]
 [-3.40e-02  3.00e-03  1.00e-03  0.00e+00]
 [-2.50e-02  4.00e-03  0.00e+00  0.00e+00]]
(4, 4)
159

　以上で、OLSによるパラメータ推定を実装できました。

　この記事では、GWRモデルの最小二乗法をプログラムで確認しました。次の記事では、図で確認します。