非線形最小ニ乗法　ニュートン法、（2乗法、二乗法）

非線形最小ニ乗法

線形の最小二乗法は、検索すると沢山ヒットしますが、
非線形については、なかなかヒットしない。
そこで、非線形の最小二乗法による、一例を作ってみました。
方法的には、ニュートン方になるのかな、エクセルのソルバーを使えば、非線形も、解けるので、
単に、どんな方法で、求められるか、興味ない人は、他のページに移った方が良いです。

		まず、解く前の基本式です。近似する方程式を左記のようにします。今回は、求める最小二乗法の式の未知のパラメータをa,b,cとします。ｘは、データです。
		最小二乗法の基礎式を左記の様にします。
		データyと、求める式の計算値ｆの差をｒとすると最小二乗法は、ｒのニ乗の∑が最小となる、a,b,cを求めることです。

そこで、ニュートン法で、求めようとすると、
テイラー展開を、利用して、いろいろな式を誘導して、その結果、
以上の行列を解いて、Δa,Δb,Δcを求めて、真の値に近づけていけば良いことになります。

		行列の計算は、エクセルのVBAでは、ワークシート関数で、逆行列と行列の積がサポートしているので非常にシンプルなプログラムになります。

実際に例を挙げます

C	q
1.6	171
4.52	228
6.8	258
8.16	284
11.5	321
12.7	335
18.2	378
29	401
38.9	410
57.3	429

左記のデータを以下の式で近似します
また基本データ　Cをx,ｑをｙとします

この式のa,b,cを、非線形最小二乗法で求めます。
ここで、雑感なのですが、いろいろvbaによる数値計算のページがあるのでが、
何故か、逆行列や行列の積のワークシート関数を利用している例は、少ないです。
自分で、ガウスの方法等で、連立方程式を解いているのですよね。
エクセルのワークシート関数を使えば、
プログラムも、シンプルになるのですが。
何ゆえに嫌われているのか、不思議です

Option Explicit


Public Function drda(x As Double, a As Double, b As Double, c As Double) As Double

drda = b * x ^ (1 + c) / (1 + a * x ^ c) ^ 2

End Function

Public Function drdb(x As Double, a As Double, b As Double, c As Double) As Double

drdb = -x / (1 + a * x ^ c)

End Function

Public Function drdc(x As Double, a As Double, b As Double, c As Double) As Double

drdc = b * a * x ^ (1 + c) * Log(x) / (1 + a * x ^ c) ^ 2

End Function

Public Function f(x As Double, a As Double, b As Double, c As Double)

f = b * x / (1 + a * x ^ c)

End Function

Public Function r(ff As Double, x As Double, a As Double, b As Double, c As Double)

r = ff - f(x, a, b, c)


End Function

Public Sub newton(a As Double, b As Double, c As Double)

Dim aa(1 To 3, 1 To 3) As Double
Dim x As Variant
Dim y(1 To 3, 1 To 1) As Double
Dim xd As Variant, yd As Variant

Dim mi As Integer
Dim i As Integer, j As Integer
Dim dn As Integer


xd = Array(1.6, 4.52, 6.8, 8.16, 11.5, 12.7, 18.2, 29, 38.9, 57.3)
yd = Array(171, 228, 258, 284, 321, 335, 378, 401, 410, 429)


For mi = 1 To 100

For i = 1 To 3
 For j = 1 To 3
    aa(i, j) = 0
 Next
 y(i, 1) = 0
 Next


For dn = 1 To 10

aa(1, 1) = aa(1, 1) + drda(CDbl(xd(dn - 1)), a, b, c) * drda(CDbl(xd(dn - 1)), a, b, c)
aa(2, 1) = aa(2, 1) + drda(CDbl(xd(dn - 1)), a, b, c) * drdb(CDbl(xd(dn - 1)), a, b, c)
aa(3, 1) = aa(3, 1) + drda(CDbl(xd(dn - 1)), a, b, c) * drdc(CDbl(xd(dn - 1)), a, b, c)

aa(1, 2) = aa(1, 2) + drdb(CDbl(xd(dn - 1)), a, b, c) * drda(CDbl(xd(dn - 1)), a, b, c)
aa(2, 2) = aa(2, 2) + drdb(CDbl(xd(dn - 1)), a, b, c) * drdb(CDbl(xd(dn - 1)), a, b, c)
aa(3, 2) = aa(3, 2) + drdb(CDbl(xd(dn - 1)), a, b, c) * drdc(CDbl(xd(dn - 1)), a, b, c)

aa(1, 3) = aa(1, 3) + drdc(CDbl(xd(dn - 1)), a, b, c) * drda(CDbl(xd(dn - 1)), a, b, c)
aa(2, 3) = aa(2, 3) + drdc(CDbl(xd(dn - 1)), a, b, c) * drdb(CDbl(xd(dn - 1)), a, b, c)
aa(3, 3) = aa(3, 3) + drdc(CDbl(xd(dn - 1)), a, b, c) * drdc(CDbl(xd(dn - 1)), a, b, c)

y(1, 1) = y(1, 1) + r(CDbl(yd(dn - 1)), CDbl(xd(dn - 1)), a, b, c) * drda(CDbl(xd(dn - 1)), a, b, c)
y(2, 1) = y(2, 1) + r(CDbl(yd(dn - 1)), CDbl(xd(dn - 1)), a, b, c) * drdb(CDbl(xd(dn - 1)), a, b, c)
y(3, 1) = y(3, 1) + r(CDbl(yd(dn - 1)), CDbl(xd(dn - 1)), a, b, c) * drdc(CDbl(xd(dn - 1)), a, b, c)

Next


x = WorksheetFunction.MMult(WorksheetFunction.MInverse(aa), y)


a = a - x(1, 1)
b = b - x(2, 1)
c = c - x(3, 1)

Debug.Print mi, a, b, c

Next



End Sub

Public Sub matrix()

Dim amin As Double, bmin As Double, cmin As Double
Dim ia As Double, ib As Double, ic As Double
Dim n As Integer
Dim ymin As Double
Dim y As Double


amin = 0
bmin = 150
cmin = 0

ymin = yy(amin, bmin, cmin)

For ia = 0 To 2 Step 0.1
 For ib = 200 To 240
  For ic = 0 To 2 Step 0.1
  y = yy(ia, ib, ic)
  If y < ymin Then
  ymin = y
  amin = ia
  bmin = ib
  cmin = ic
  
  Debug.Print amin, bmin, cmin, y
  
  End If
  Next
 Next
Next

Debug.Print amin, bmin, cmin

Call newton(amin, bmin, cmin)

End Sub

Public Function yy(a As Double, b As Double, c As Double) As Double
Dim xd As Variant, yd As Variant
Dim y As Double
Dim i As Integer

y = 0

xd = Array(1.6, 4.52, 6.8, 8.16, 11.5, 12.7, 18.2, 29, 38.9, 57.3)
yd = Array(171, 228, 258, 284, 321, 335, 378, 401, 410, 429)

For i = 1 To 10

y = y + r(CDbl(yd(i - 1)), CDbl(xd(i - 1)), a, b, c) ^ 2


Next

yy = y


End Function

まず関数ｒのa,b,cで偏微分した式を
定義しました。

求める関数fを定義しました

最小二乗法に使うrの式を定義しました。

この部分は、ニュートン法の関数の開始点です。
初期値をa,b,cを引数にしています。

行列に値の代入しています。

逆行列を算出して、Δa,Δb,Δcを求めるために行列の
積を計算しています

Δa,Δb,Δcの値が、xに現るので、
初期値のa,b,cを補正して、
新たな、初期値a,b,c,を求めました

一度の計算による、近似値の値の出力
収束しない場合等を考えて、
近似計算を100繰り返しています。
収束したら、プログラムをストップしたかったら
条件文等で、止めてもいいと思います。

今回の初期値を求めるために
a,b,cの取りうる値を推定して、
∑r^2が、最小値になる
求める解の近傍のa,b,cを求めています。

∑ｒ＾２の関数です

上記プログラムコードをコピーペーストで、貼り付けたら以下のコマンドで、実行しました。

計算結果は以下の様になりました。

10回程度の計算で収束していることがわかります。

次に、関数ｆを解析的に偏微分するのが難しかった時に、数値微分で、偏微分した場合の修正部分と
結果を表します。

Const ca As Double = 0.001
Const cb As Double = 0.01
Const cc As Double = 0.001

Public Function drda(x As Double, a As Double, b As Double, c As Double) As Double

drda = 1# / (12# * ca) * (-f(x, a - 2# * ca, b, c) + 8# * f(x, a - ca, b, c) - 8# * f(x, a + ca, b, c) + f(x, a + 2# * ca, b, c))

End Function

Public Function drdb(x As Double, a As Double, b As Double, c As Double) As Double

drdb = 1# / (12# * cb) * (-f(x, a, b - 2# * cb, c) + 8# * f(x, a, b - cb, c) - 8# * f(x, a, b + cb, c) + f(x, a, b + 2# * cb, c))

End Function

Public Function drdc(x As Double, a As Double, b As Double, c As Double) As Double

drdc = 1# / (12# * cc) * (-f(x, a, b, c - 2# * cc) + 8# * f(x, a, b, c - cc) - 8# * f(x, a, b, c + cc) + f(x, a, b, c + 2# * cc))

End Function

数値微分は、5点法で行ないました。　計算結果は以下のようになりました。

極端な、精度を求めなければ数値微分で、
計算しても、問題ないような感じです。
でも、数値微分のΔxをうまく調整しないと
駄目なようです。

気が付く人もいるかと思いますが。
これって、線形の最小二乗法にも使えないかという、疑問が生まれます。
そうなのです。非線形の方法でも、線形の方法としても、使えるのです。

最小ニ乗法　最小2乗法　最小二乗法