◆ Covar()

Sub Covar	(	N As	Long,
		R() As	Double,
		Ipvt() As	Long,
		Tol As	Double,
		Info As	Long
	)

非線形最小二乗近似の分散共分散行列

目的: 本ルーチンは, Lmder, Lmder1, Lmdif, LmstrまたはLmstr1により求められた非線形最小二乗近似の分散共分散行列を求める.

M x N 行列 A について, 分散共分散行列と呼ばれる N x N 対称行列 C が次のように定義される.
C = (A^T*A)^(-1)

本ルーチンは, Lmder, Lmder1, Lmdif, LmstrまたはLmstr1により得られるQR分解結果を使ってこれを求める.

引数

[in]	N	分散共分散行列の次数. (N > 0)
[in,out]	R()	配列 R(LR1 - 1, LR2 - 1) (LR1 >= N, LR2 >= N) [in] QR分解 AP = QR より得られる上三角行列 R. (P は Ipvt()で与えられる置換行列, Q は M x N 直交行列である) [out] 分散共分散行列 (対称な正方行列).
[in]	Ipvt()	配列 Ipvt(LIpvt - 1) (LIpvt >= N) 置換行列を定義するピボット情報.
[in]	Tol	行列 A のランク数を決めるための許容値. (Tol > 0) R の対角成分のうち, (絶対値最大の対角成分×Tol)より絶対値が小さい対角成分があればランク落ちと見なす.
[out]	Info	= 0: 正常終了. = -1: パラメータ N の誤り. (N <= 0) = -2: パラメータ R() の誤り. (配列R()の大きさが不足) = -3: パラメータ Ipvt() の誤り. (配列Ipvt()の大きさが不足) = -4: パラメータ Tol の誤り. (Tol <= 0)

出典: netlib/minpack

使用例: 次のデータをモデル関数 f(x) = c1*(1 - exp(-c2*x)) で近似する.
2つのパラメータc1, c2をLmder1により定め, Covarによりその分散共分散行列を求める.
f(x) x

10.07 77.6

29.61 239.9

50.76 434.8

81.78 760.0

Sub FLmder(M As Long, N As Long, X() As Double, Fvec() As Double, Fjac() As Double, IFlag As Long)

Dim Xdata(3) As Double, Ydata(3) As Double, I As Long

Ydata(0) = 10.07: Xdata(0) = 77.6

Ydata(1) = 29.61: Xdata(1) = 239.9

Ydata(2) = 50.76: Xdata(2) = 434.8

Ydata(3) = 81.78: Xdata(3) = 760

If IFlag = 1 Then

For I = 0 To M - 1

Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))

Next

ElseIf IFlag = 2 Then

For I = 0 To M - 1

Fjac(I, 0) = Exp(-Xdata(I) * X(1)) - 1

Fjac(I, 1) = -Xdata(I) * X(0) * Exp(-X(1) * Xdata(I))

Next

End If

End Sub

Sub Ex_Covar()

Const M = 4, N = 2

Dim X(N - 1) As Double, Fvec(M - 1) As Double, R(M - 1, N - 1) As Double

Dim Tol As Double, Ipvt(N - 1) As Long, Info As Long

Dim S As Double, I As Long

Tol = 0.00000001 '1.0e-8

X(0) = 500: X(1) = 0.0001

Call Lmder1(AddressOf FLmder, M, N, X(), Fvec(), R(), Tol, Ipvt(), Info)

Debug.Print "Info =", Info

Call Covar(N, R(), Ipvt(), Tol, Info)

For I = 0 To M - 1

S = S + Fvec(I) ^ 2

Next

S = S / (M - N)

Debug.Print "Cov ="

Debug.Print R(0, 0) * S, R(0, 1) * S

Debug.Print R(1, 0) * S, R(1, 1) * S

Debug.Print "Info =", Info

End Sub

Covar
Sub Covar(N As Long, R() As Double, Ipvt() As Long, Tol As Double, Info As Long)
非線形最小二乗近似の分散共分散行列

Lmder1
Sub Lmder1(F As LongPtr, M As Long, N As Long, X() As Double, Fvec() As Double, Fjac() As Double, Tol As Double, Ipvt() As Long, Info As Long, Optional Info2 As Long)
非線形最小二乗法 (レーベンバーグ・マルカート法) (シンプルドライバ)

実行結果: Info = 0

Cov =

20.5868644246757 -5.52380489748924E-05

-5.52380489748924E-05 1.48678858325168E-10

Info = 0