誤差論と最小二乗法

第8回 線形モデル – その2

前回は基本の線形モデルを取り上げ、最小二乗法により未知パラメータを推定しました。今回は、もう一つの推定法である最尤法で未知パラメータを求めてみます。次に線形モデルを少し一般化したものを考えます。具体的には、image1の場合です。各観測値の分散が異なる場合や相関を持つ場合に相当します。

 

1. 最尤(推定)法

最小二乗法は線形モデルにおける基本的な方法ですが、様々なデータの中には線形モデルで記述するには不適切なものもあります。そのような場合に一般に用いられる方法が最尤法です。

 

最尤法は第6回に紹介しました。尤度関数を最大にする推定値を求める方法です。尤度関数とは観測値の確率密度関数において観測値は観測された値に固定し、パラメータを変数とみなしたものです。ここでは最尤法を正規分布の線形モデルに用いてみましょう。

 

線形モデル

image2

において、観測値の分布を(n次元の)正規分布(付録B参照

 

image3

と仮定すると(image4に注意)、未知パラメータは、尤度関数

 

image5

あるいは、対数尤度

 



を最大にするもの(最尤推定値)として求まります。

 

image7の最尤推定値は上式の第3項を最大にするものですから残差二乗和を最小にするもの、つまり最小二乗法と同じになります。

 

最尤法は様々な分布を持つデータに用いることができる一般的な方法ですが、定義からもわかるように観測値の分布関数を仮定する必要があります。最小二乗法では特に分布関数の形に仮定はなく、分散(あるいは重み)のみが仮定されます。

 

2. 一般の最小二乗法

線形モデルにおいて各観測値の分散が必ずしも同じでない場合を考えます。すると、image9を測定データ、image10を既知の計画行列、image7を未知パラメータ、image11をランダム誤差とするとモデルは

 

image12

 

と書けます。ここで、image13は既知の正定値行列、image14は一般には未知とします。式image15は、観測値を理論式と誤差で表すので観測方程式と呼ばれます。

 

このとき付録A2より

image16

 

と分解できるので、image9image17を左から掛けると(2)は、

 

image18

 

image19

 

と変形され基本モデルと同じ形になります。

 

最小二乗の条件は、

 

が最小なることと言い換えられます。ここで

 

image22

 

は、重み行列と呼ばれます。従って、最小二乗の条件は単なる二乗和=長さの二乗ではなく、重みが付いた長さの二乗=重み付き残差二乗和:image23が最小になることです。

正規方程式は、前回の式(8)においてimage9image24image10image25に置き換えて

 

image26

 

となり、最小二乗解は、

 

 

と求まります。この解は普通の最小二乗解と同じくBLUEとなります。

 

また、image28の誤差行列とimage14の推定値も同様に

 

image29

 

となります。

 

重みと分散
ある観測image31の重みimage32は分散image33に逆比例する量として定義されます。

 

image34

 

分散が小さければその観測の重みは大きくなります。特に、各観測値image35がおのおの独立で異なる分散image33を持つとすると

 

image36

 

より、重み行列は

 

image37

 

となります。image38としたものが(2)で、基準となる観測image31を選びimage39とするとimage40となりますからimage14は単位重量当たりの分散あるいは基準分散と呼ばれます。時には各観測値の分散が既知の場合があります。その場合は、image41として、image42とすればよいのですが、image14の値は最小二乗解には影響を与えません。

 

3.いくつかの例

ここで測量に関連した例を挙げてみたいと思います。

 

例1.水準測量
簡単な水準路線で未知点の標高を求めてみましょう(図1)。既知点image44の標高がimage45、未知点image46image47の標高が未知パラメータimage48、測定値は比高image49です。

 

image43図1.水準路線

 

 

モデルは、

image50

image51

image52


ですが、これを

 

image53

 

の形に行列で書くと、

 

image54

 

として、

 

image55

 

となります。

水準測量の場合、観測の重みは路線長に反比例するという経験則があるのでそれを採用すると、重量行列は各路線長をimage56として

 

image57

です。(5)、(6)の各行列を計算すると

 

image58

image59

から(image60の逆行列の計算は付録A2参照)、最小二乗解が

 

image61

 

と求まります。

 

この例は簡単な水準網でしたが、どのような水準網も観測値(比高)は未知パラメータである水準点標高の和や差で表せるので線形モデルになっています。

 

例2.基準点測量(距離測量による平面位置決定)

基準点測量は、距離、角度、GNSSなどの観測を使って新しい点の位置を求める測量です。点の座標と観測値は一般に非線形の関係ですから、テイラー展開(付録C参照)で線形化して(2)の形にすることで最小二乗法を適用します。

 

簡単な距離測量による平面位置の決定の問題(図2)について見てみましょう。

 

image62

図2.距離観測による位置決定

 

 

image63が既知点、image64が未知点でimage65からimage66までの距離image67を測定したとします。観測の重みはそれぞれimage65とします。

 

距離は位置座標によって次のように非線形関数で表されます。

 

image68

 

未知座標を概算値(添字0)+補正量: image69としてテイラー展開すると(微分の詳細は略)、

 

image70

 

となるので、各距離の観測方程式は、

 

image71

行列でまとめて書くと、

image72

 

image73

 

image74

 

となります。最小二乗解は例1と同様に計算できPの位置ベクトルが

 

image75

と求まりますが、注意点が一つあります。観測方程式はあくまで非線形関数を近似した式のため、最初に解いた解には省略した項の影響が残っている可能性があります。従って、最初の解を概略値としてもう一度最小二乗解を求め、さらにこれを何回か繰り返し収束したものを最終的な解としています。測地網の場合、数学モデルが堅固なので普通、1回の反復改良計算で収束します。ただし、初期値と推定値が大きい場合、2回程度の反復改良計算が必要なことがあります。

 

例3.GNSS測量(二重位相差による相対測位)

受信機を点A(既知点) とB(未知点)におき、Bの位置を位相観測による相対測位で決定する測量を考えましょう。精密測量では一般に二重位相差をとって時計誤差を消去します(図3)。

 

image76図3.二重位相差の観測

 

 

ここではABを短基線とし、大気による誤差も無視できるとすると単純化された観測方程式は次のようになります。

 

image77

 

ここで、image78は二重位相差、

 

image79

 

image80は搬送波波長、

 

image81は整数アンビギュイティー

 

です。実際は、このような観測が衛星のペア×エポック数だけあります。 image82は受信機と衛星間の距離で位置座標の線形関数ではありませんから、例2と同様に線形化します。例えば、4衛星1エポックの観測からは3つの独立な二重位相差が得られ(図4)、未知数はBの座標(の補正量)と3つのアンビギュイティーの6つです:image83

 

image84

図4.4衛星1エポックの二重位相差

 

 

image85を線形化すると、計画行列image10image86番目のエポックに対応する行は、次のようになります。

 

image87

 

微分を行うと、

 

image88

image89

 

などから、例えば

 

image90

 

となります。行の数は観測の数と同じです。このようにして線形化された式を最小二乗法で解くわけですが、例1,2の場合と違うことがあります。それは、観測データ(二重位相差)の間に相関があることです。図4の二重位相差の式を見ると、一重位相差image91が3つに共通に含まれていることからもわかります。各点での位相観測は独立と考えられ共分散行列は対角行列ですが、二重位相差は位相の線型結合になっています。二重位相差の観測値ベクトルimage92は、位相観測値ベクトルimage93を行列image94で変換して,

 

image95

 

と得られる(image94の要素は+1、1、0のいずれかです)ので、image92の共分散行列は誤差伝播(前回参照)より、

 

image96

となります。従って重量行列は、

 

image97

 

と書け、この重量行列を用いて最小二乗解を求めることになります。

 

 『第8回付録 線形(線型)代数の基礎2』へ

Share on FacebookTweet about this on TwitterShare on LinkedInEmail this to someonePrint this page