はじめに
『Pythonで学ぶ空間データサイエンス入門』の独学ノートです。本の内容から寄り道してアレコレ考えます。
本を読んだ上で補助的に読んでください。
この記事では、GWRモデルの最小二乗法について、数式を使って解説します。
【前の内容】
【他の内容】
【今回の内容】
4.3.4-7 GWRモデルの最小二乗法の導出
地理的加重回帰モデル(GWRモデル・geographically weighted regression model)に対する最小二乗法(OLS・ordinary least square)におけるパラメータの計算式を導出します。
線形回帰モデルや最小二乗法については「3.2.1-3:線形回帰モデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。
モデルの設定
まずは、GWRモデルの定義を数式で確認します。
個の地点について、各地点 の説明変数(入力データ)を 次元ベクトル とします。
個の説明変数をまとめて、 の行列 とします。
各変数に関して、 番(列)目の要素を とします。(インデックスを0から割り当てているので) 次元(列)です。
番目の要素 に対する基準地点 の真の係数パラメータを として、 個の真の係数パラメータをまとめて、 次元ベクトル とします。
番目の要素 は定数項に対応します。
各地点 の誤差項を として、 個の誤差項をまとめて、 次元ベクトル とします。
各誤差 は、平均が で分散が (標準偏差が )の正規分布に従って、独立に生成されると仮定します。
番目の誤差項 の共分散を とすると、誤差項間は無相関(無関係に生成される)なので、 です。各誤差項自身との共分散は分散 です。
基準地点 の重み行列を の行列 とします。
異なる地点 の組み合わせに対応する要素(非対角要素)は です。よって、重み行列は対角行列です。
各地点(同じ地点の組み合わせ) に対応する要素(対角要素)は で表され、重み関数(カーネル)に従い、基準の地点 との距離に応じて値が決まります。基準の地点 に対応する要素は です。
GWRモデルの重み行列については「【Python】4.3.3:GWRモデルの重み行列の可視化【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。
各地点 の被説明変数(出力データ)を として、 個の被説明変数をまとめて、 次元ベクトル とします。
(1つの地点における)基準地点を とするGWRモデルは、次の式で定義されます。
途中式の途中式(クリックで展開)
- 1: 線形回帰モデルの式 の両辺に重み を掛けます。
- 2: 括弧を展開します。
- 3: ベクトルの要素を明示します。
- 4: ベクトルの内積を計算します。
- 5: 個の項の和を でまとめます。
- 6: (項の形が異なる) 番目の項 を取り出し( 個の要素の和を でまとめ)ます。
個のデータをまとめると、次の式で計算できます。
途中式の途中式(クリックで展開)
- 1: 線形回帰モデルの式 の両辺に左から重み を掛けます。
- 2: 括弧を展開します。
- 3: 行列とベクトルの要素を明示します。
- 4-5: 行列とベクトルの積や行列の積を計算します。
- 6: ベクトルの和を計算します。
ルートをとった(2分の1乗の)重み を用いることで、後ほど2乗される際に式(表記)がシンプルになります。数値上は、ルートをとった重み関数 により重みを設定して、1乗の重み を用いてモデルを定義するのと同じ操作です。
ここまでは、GWRモデルを確認しました。次からは、位置(距離)データから設定される重み行列 と観測データ が得られたときのGWRモデルのパラメータを推定します。(真のパラメータが基準地点に依存するなら被説明変数が基準地点ごとに変わってしまうのでは?定数項で調整されることになるのかな?)
係数パラメータの計算式
GWRモデルに対するOLSによる係数パラメータの推定値の計算式を導出します。
説明変数の 番目の要素 に対する基準地点 の係数パラメータの推定値(最小二乗推定量)を として、 個の係数パラメータの推定値をまとめて、 次元ベクトル とします。
基準地点を とするときの各地点 の被説明変数の推定値を として、 個の被説明変数の推定値をまとめて、 次元ベクトル とします。
各被説明変数 の推定値 は、推定パラメータ を用いて求めた値です。
途中式の途中式(クリックで展開)
- 1: 線形回帰モデルの式 の両辺に重み を掛けます。
- 2: ベクトルの要素を明示します。
- 3: ベクトルの内積を計算します。
- 4: 個の項の和を でまとめます。
- 5: (項の形が異なる) 番目の項 を取り出し( 個の要素の和を でまとめ)ます。
個の重み付けした被説明変数の推定値をまとめると、次の式で計算できます。
途中式の途中式(クリックで展開)
- 1: 線形回帰モデルの式 の両辺に左から重み を掛けます。
- 2: 行列とベクトルの要素を明示します。
- 3-4: 行列とベクトルの積や行列の積を計算します。
観測データとして得られた を実現値、推定パラメータから求めた を理論値とも呼びます。
基準地点を とするときの各地点 に関して、重み付けした実現値 と理論値 の差を残差 とします。
個の(1つの基準地点 による)各地点の残差をまとめると、次の式で計算でき、 次元ベクトル とします。
途中式の途中式(クリックで展開)
- 1: 重み付けした実現値 の式(1)と理論値 の式(2)の差の式を立てます。
- 2: を括り出します。
- 3: 行列とベクトルの要素を明示します。
- 4: ベクトルの差を計算します。
- 5: 行列とベクトルの積を計算します。
残差 の式に関して、理論値 の式(1)で置き換えて使います。
残差 の2乗和(残差平方和)を考えます。
途中式の途中式(クリックで展開)
- 1: の式(3)より、 の内積の式を立てます。
- 2: 括弧を展開します。
- 3: 転置行列の性質 より、括弧を展開します。
- 3-4: 重み行列は対角行列なので、、また です。
- 4: スカラは転置しても変化しないので、2つ目の項は で置き換えられ、転置行列の性質 より、 となります。
- 5: 式を整理します。
残差平方和 を推定パラメータ に関して微分します。
途中式の途中式(クリックで展開)
- 1: に対する の偏微分の式を立てます。 に関する微分なので、 は定数として扱います。
- 2: 和の微分 より、項ごとの微分の和に分割します。
- 3: 定数の微分 より、1つ目の項は0ベクトルになり消えます。
- 3: 定数倍の微分 より、 の係数を の外に出します。
- 3: 対象行列の二次形式の微分 より、 を1つの行列とみなして変形します。
- 4: 変数の微分 より、2つ目の項の微分が単位行列になり消えます。
残差平方和の微分 が0ベクトル となる推定パラメータ を求めます。
途中式の途中式(クリックで展開)
- 1: とおきます。
- 2-3: について式を整理します。
- 3: 両辺に左から の逆行列を掛けます。
の計算式が得られました。
以上で、GWRモデルの係数パラメータの推定値の計算式が求まりました。
分散パラメータの計算式
GWRモデルに対するOLSによる誤差項の生成分布の分散パラメータの推定値の計算式を導出します。
各基準地点 に関して、地点 の重み付けしていない残差 を考えます。
途中式の途中式(クリックで展開)
- 1: 実現値 と理論値 の差の式を立てます。
- 2: 線形回帰モデルの式 で置き換えます。
- 3: の式(4)で置き換えます。
- 4: 線形回帰モデルの式 、 で置き換えます。
- 5: 括弧を展開します。
- 6: 逆行列の性質 より、 が単位行列になり消えます。
- 7: 式を整理します。
後の項の横ベクトルについて、 の転置とおきます。
逆行列部分は、、、 の行列の積なので、 の行列です。 は、、、、 の行列の積なので、 のベクトルになります。
残差 の式(6)を で置き換えます。
基準地点 のときの地点 の残差(観測値と推定値の差) と誤差(観測値と真値の差) の関係式が得られました。
個の各基準地点の残差をまとめると、次の式で計算でき、 次元ベクトル とします。
途中式の途中式(クリックで展開)
- 1: 個の基準地点に関して、 の式(7)を縦方向に結合した式を立てます。
- 2: ベクトルの差の計算を分割します。
- 3: 行列とベクトルの積の計算を分割します。
後の項の行列について、 とおきます。
は、 個の 次元ベクトル を行方向に結合した行列なので、 の行列になります。
残差 の式(9)を で置き換えます。
各基準地点の残差(観測値と推定値の差) と誤差(観測値と真値の差) の関係式が得られました。
同様に、 の式(5)を用いた を考えます。
途中式の途中式(クリックで展開)
- 1: 個の基準地点に関して、 の式(5)を の式(6)で置き換えて、縦方向に結合した式を立てます。
- 2: ベクトルの差の計算を分割します。
- 3: 行列とベクトルの積の計算を分割します。
- 4: の式(10)で置き換えます。
- 5: を括り出します。
( は観測できる値であるが は観測できない値なので、) の計算式が得られました。
残差 の式(11)を用いて、残差平方和を考えます。
途中式の途中式(クリックで展開)
- 1: の式(11)より、 の内積の式を立てます。
- 2: 転置行列の性質 より、 となります。
- 3: 括弧を展開します。
- 4: 式を整理します。
- 5: 括弧を展開します。
さらに、要素を明示して計算します。
途中式の途中式(クリックで展開)
1つ目の項は
となります。
2つ目の項は
となります。
3つ目の項について
とおくと
となります。
残差平方和 の期待値を考えます。
途中式の途中式(クリックで展開)
- 1: の期待値の式を立てます( の式(12)の両辺で期待値をとります)。
- 2: 期待値の性質 より、和の期待値を期待値の和に分割します。
- 2: 期待値の性質 より、定数倍の期待値を期待値の定数倍に変形します。
- 3: によって定義される は と無関係なので、 の各要素 を の外に出せます。
各誤差項 は平均 、分散 の正規分布に従い独立に生成されると仮定しているので、 の期待値は であり、 の積の期待値は の共分散 となります。また、同じ組み合わせ の共分散は の分散 です。
よって、次の関係が得られます。
詳しくは「線形回帰モデルの最小二乗法の導出」を参照してください。
この関係を用いて、残差平方和の期待値 の式を整理します。
途中式の途中式(クリックで展開)
- 1-2: 各組み合わせの積の期待値 を式(13)で置き換えます。
- 2: 3つの項それぞれで、 に関する項の和 を展開すると に関する項のみが残ります。
- 3: 定数の総和は定数倍 です。
- 3: は各データと無関係なので の外に出せます。
- 4: を括り出します。
- 5: は 、 は の対角要素の和なので、行列のトレース(対角要素の和) に置き換えます。
この式を分散パラメータ について解きます。
の計算式が得られました。
残差平方和の期待値 を残差平方和 で近似
して、置き換えます(なぜこんな操作をしてもいいのか?)。
ここで、次のようにおきました。
を有効パラメータ数と言います。
以上で、線形回帰モデルの誤差項の分散パラメータの計算式が求まりました。
係数パラメータの分散共分散行列の計算式
GWRモデルに対するOLSによる係数パラメータの推定値の分散共分散行列の計算式を導出します。
誤差項 の期待値は、 より、0ベクトルになるのでした。
また、誤差項 の分散共分散行列は、 の式(13)より、対角要素が分散パラメータ の対角行列になるのでした。
詳しくは「線形回帰モデルの最小二乗法の導出」を参照してください。
係数パラメータの推定値 の式(4)を変形します。
途中式の途中式(クリックで展開)
- 2: 線形回帰モデルの式 で置き換えます。
- 4: 括弧を展開します。
- 5: 逆行列の性質 より、 が単位行列になり消えます。
推定パラメータ の期待値は、真のパラメータ になります。
途中式の途中式(クリックで展開)
- 1: の期待値の式を立てます( の式(16)の両辺で期待値をとります)。
- 2: 和の期待値 より、項ごとの期待値に分割します。
- 3: 観測値 は定数として扱うので、 の外に出せます。
- 4-5: の式(14)で置き換えると、後の項が0ベクトルになり消えます。
よって、推定パラメータ と推定パラメータの期待値 の差は、次の式で表わせます。
途中式の途中式(クリックで展開)
- 1: の差の式を立てて、 の式(17)で置き換えます。
- 2: の式(16)で置き換えます。
- 3: 式を整理します。
推定パラメータ の分散共分散行列を考えます。
途中式の途中式(クリックで展開)
- 1: 分散共分散行列の定義 より、 の分散共分散行列の式を立てます。
- 1-2: の差を式(14)で置き換えます。
- 3: 転置行列の性質 、転置行列と逆行列の関係 、重み行列は対角行列なので より、 となります。
- 4: 観測値 は定数として扱うので、 の外に出せます。
- 5: を式(15)で置き換えます。
- 6: スカラ を前に出します。
この式の行列について、 とおきます。
逆行列部分は の行列でした。 は、、、 の行列の積なので、 の行列になります。
ちなみに、、 の関係です。
推定パラメータ の分散共分散行列の式(19)を で置き換えます。
の分散共分散行列の計算式が得られました。
この記事では、GWRモデルの最小二乗法を数式で確認しました。次の記事では、プログラムで確認します。
参考文献
おわりに
最後惜しいなぁ、真ん中の重み行列の1つを消すなり移動するなりできれば、線形回帰のときのようにスッキリした式になるのになぁ。気になるなぁ。
ところでハット行列とはなんですか。ここでハット行列と呼ばれているものと3章で導出したものとは別ですよね。謎です。
基準地点 に対してなのか(基準地点も含めた)ある地点 に対してなのか、はたまた 個の各地点 に対してなのかを意識して読む(変形する)必要があるのですが、途中から 個の各基準地点 に対して操作(計算)していて、そこに気付くまで随分苦労しました。
異なる基準地点の重みやパラメータを扱っている場合に区別するため(本にはない記号を導入して) のように基準地点 を明示的に表記しました。ベクトルの要素らしさ(行列の要素らしくなさ)のため のように基準地点番号と地点番号を分けて表記することも考えたのですが、不必要に記号が煩い気がして止めました(他のシリーズではたまに使っています)。自分で書いているとどう表記しても分かりやすく思ってしまうので悩みます。
重みやパラメータ についても、(基準地点 に依存して決まるのを関数的に表しているのだと思うのですが、)煩いので のように簡素化して表記しました。
(この記事では登場させずに済んだのですが)本や他の記事では、 個の基準地点に関してまとめて、行列 や として扱う(と行列計算により表記や処理が楽になる)こともあるので、基準地点番号を行番号として用いて行列 や のように流用できる表記にしました。
とすると重みは のようにした方が行列の成分が分かりやすいと思いますが、他との兼ね合いで基準地点番号を1つ目の添字で示すということで納得することにしました。
基準地点なのか各地点なのかや、どの要素がどの要素に影響するのかを明示したくて、ベクトルや行列の計算を要素レベルの計算まで何度か行いました。清書し終わってから気付いたのですが、次の式変形ではベクトルや行列の状態の計算に戻ってから行っており、三歩進んで二歩下がる的な展開が繰り返されいて分かりにくい書き方になっていますね。
これも自分で書いていると意図を知っているので分かりやすく思ってしまうがゆえです。
ただGWRでは、同じものを示しつつ異なる意味を持つ添字(地点番号)が各所に出てくるので、要素レベルでの対応関係を確認しておくと理解しやすいと思います。読み取りにくいけどね。
そんなこんなでなんとか理解できたと思われます。
分かってくると次の疑問が生じて、今回の内容も最小二乗法なわけですよね、線形回帰モデルの最小二乗法をOLSと呼んで地理的加重回帰モデルの最小二乗法をGWRと呼ぶのが気になってしょうがないです。そもそも今回の操作を重み付き(残差の)最小二乗法と呼びたい気もします。
2024年9月29日は、モーニング娘。10期メンバーの加入13周年の日です。
あゆみんの卒コンの日程も決まりカウントダウンがより明確に聞こえてきますが(泣)、あゆみんがラストツアー共々晴れやかな気持ちで楽しめますように。
【次の内容】