Avenueで最小二乗法 (2003.10.15)
技術資料トップへもどる
理論

測定値yを説明する式f(x)があり、n個の式

・・・(1)
が得られているものとします。ここで測定値yには誤差vを含むものとして、測定値から理論値を差し引いた残差の自乗和牌2を最小にするという条件のもとに、未知の係数a(k), a(k-1),・・・a(1), a(0)を最小二乗法的に求めます。

式(1)のyをy+vでおきかえ、yを右辺に移して行列表現で書き換えると、(1)は

・・・(2)
ここに、

となります(式の導出の詳細は参考書などを参照してください)。

tV・V=最小(tVVの転置行列)の条件から、正規方程式

・・・(3)
が得られます。ここでBに逆行列が存在すれば (3)は解けて、

・・・(4)
が得られます。


参考書

測量計算法−行列最小二乗法から網平均まで、原田健久、鹿島出版会(2001)


数値例

理論式を y = ax + b として、6回の測定値から未知の係数a、bを最小二乗法で求めます。

y x
1.40.9
1.81.6
2.51.9
2.12.5
2.42.9
3.03.4
6対の測定値で式を書き下すと、
1.4 = 0.9a + b
1.8 = 1.6a + b
2.5 = 1.9a + b
2.1 = 2.5a + b
2.4 = 2.9a + b
3.0 = 3.4a + b
となります。これを行列表現にすると
ただし、 となります。

Avenueスクリプト

以下は行列演算ライブラリを使って上の問題を解くAvenueスクリプトです。 実行すると、



となり、a=0.543269、b=1.00481が得られます。


'-------------------------------------------------------
' 最小二乗法の行列解法
'
' (C)2003 GIS Institute, Inc.
'-------------------------------------------------------
'観測方程式
'1.4 = 0.9a + b
'1.8 = 1.6a + b
'2.5 = 1.9a + b
'2.1 = 2.5a + b
'2.4 = 2.9a + b
'3.0 = 3.4a + b

'測定値
Yn = {1.4, 1.8, 2.5, 2.1, 2.4, 3.0}
Xn = {0.9, 1.6, 1.9, 2.5, 2.9, 3.4}

'----- 操作する行列を宣言する -----
'行列A:6行2列
A = av.Run("Matrix.Declare",{6,2})
'行列U:6行1列
U = av.Run("Matrix.Declare",{6,1})

'----- AとUに要素を代入する -----
'リストの添字は0から始まることに注意
for each idx in 1..Xn.Count
  'Aのidx行1列目
  av.Run("Matrix.SetElement",{A,idx,1,Xn.Get(idx-1)}) 
  'Aのidx行2列目
  av.Run("Matrix.SetElement",{A,idx,2,1})
  'Uのidx行1列目
  av.Run("Matrix.SetElement",{U,idx,1,Yn.Get(idx-1)})
end

'----- Aの転置行列を求める -----
tA = av.Run("Matrix.Transpose",A)

'----- Bを求める -----
' B = tA・A
B = av.Run("Matrix.Multiply",{tA,A})

'----- Rを求める -----
' R = tA・U
R = av.Run("Matrix.Multiply",{tA,U})

'----- Bの逆行列を求める -----
'Bは2行2列
invB = av.Run("Matrix.Inverse2",B)

'----- Xを求める -----
X = av.Run("Matrix.Multiply",{invB,R})

'----- Xの要素を参照する -----
coeffA = av.Run("Matrix.GetElement",{X,1,1})
coeffB = av.Run("Matrix.GetElement",{X,2,1})
msgbox.Info("Y = "+ coeffA.AsString + " X + " + coeffB.AsString, "最小二乗法の行列解法")

'----- 誤差の評価 -----
'理論値
theory = av.Run("Matrix.Multiply",{A,X})
'残差(観測値と理論値の差)
V = av.Run("Matrix.Subtract",{theory,U})
tV = av.Run("Matrix.Transpose",V)
'残差の2乗和
sigma = av.Run("Matrix.GetElement",{av.Run("Matrix.Multiply",{tV,V}),1,1})
'標準偏差
stdDev = ( sigma / (Xn.Count-1) ).SQRT
msgbox.Info("標準偏差: "+stdDev.AsString,"")

'----- ここまで -----

技術資料トップへもどる

©2003-2006, 有限会社GISインスティテュート. All rights reserved.
〒183-0035 東京都府中市四谷3-49-18
電話・FAX:042-405-3901
gis-i@nifty.com