◆ N2g1_r()

Sub N2g1_r	(	M As	Long,
		N As	Long,
		X() As	Double,
		Info As	Long,
		YY() As	Double,
		YYp() As	Double,
		IRev As	Long,
		Optional NFcall As	Long,
		Optional NFjcall As	Long,
		Optional Niter As	Long
	)

非線形最小二乗法 (適応アルゴリズム) (シンプルドライバ) (リバースコミュニケーション版)

目的: 本ルーチンは, ガウス・ニュートン法, レーベンバーグ・マルカート法などを組み合わせ拡張した適応アルゴリズムにより, M個のN変数非線形関数の二乗和の最小点を求める.
min Σfi(x1, x2, ..., xn)^2 (ただし, Σは i = 1 〜 M)

IRevに従って関数値およびヤコビ行列をユーザーが与える必要がある.

N2g1_rはデフォルトパラメータでN2g_rを使用するのと等価である. ただし, 分散・共分散行列と回帰診断ベクトルは計算しない.

引数

[in]	M	データ数. (M > 0)
[in]	N	パラメータ数. (0 < N <= M)
[in,out]	X()	配列 X(LX - 1) (LX >= N) [in] 初期近似解. [out] IRev = 0: 求められた解ベクトル. IRev = 1, 2: 関数値またはヤコビ行列を求める点.
[out]	Info	= 0: 正常終了. = -1: パラメータ M の誤り. (M < N) = -2: パラメータ N の誤り. (N < 1) = -3: パラメータ X() の誤り. = -5: パラメータ YY() の誤り. = -6: パラメータ YYp() の誤り. = 7: 特異収束. (近傍のヘッセ行列が特異になった) = 8: 誤収束. (誤った点での収束と思われる. 目標精度が小さすぎる可能性がある) = 9: 関数評価回数の最大値を超えた. = 10: 反復回数の最大値を超えた. = 63: X()の初期点においてF(X)を求めることができない. = 65: X()において微分係数を求めることができない.
[in]	YY()	配列 YY(LYY - 1) (LYY >= M) IRev = 1の場合, 再呼び出し時にX()における関数値を与えること.
[in]	YYp()	配列 YYp(LYYp1 - 1, LYYp2 - 1) (LYYp1 >= M, LYYp2 >= N) IRev = 2の場合, 再呼び出し時にX()におけるヤコビ行列を与えること.
[in,out]	IRev	リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: X()における関数値を求め YY()に設定すること. YY()以外の変数を変更してはならない. = 2: X()におけるヤコビ行列を求め YYp()に設定すること. YYp()以外の変数を変更してはならない.
[out]	NFcall	(省略可) 関数評価回数(IRev = 1で戻った回数).
[out]	NFjcall	(省略可) ヤコビ行列評価回数(IRev = 2で戻った回数).
[out]	Niter	(省略可) 反復回数.

出典: netlib/port

使用例: 次のデータをモデル関数 f(x) = c1*(1 - exp(-c2*x)) で近似する. 2つのパラメータc1, c2を非線形最小二乗法により定める.
f(x) x

10.07 77.6

29.61 239.9

50.76 434.8

81.78 760.0

初期値は, c1 = 500, c2 = 0.0001 とする.
Sub FN2g(M As Long, N As Long, X() As Double, Nf As Long, Fvec() As Double)

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

For I = 0 To M - 1

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

Next

End Sub

Sub JN2g(M As Long, N As Long, X() As Double, Nf As Long, Fjac() As Double)

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

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 Sub

Sub Ex_N2g1_r()

Const M = 4, N = 2

Dim X(N - 1) As Double, Info As Long

Dim YY(M - 1) As Double, YYp(M - 1, N - 1) As Double, IRev As Long

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

IRev = 0

Do

Call N2g1_r(M, N, X(), Info, YY(), YYp(), IRev)

If IRev = 1 Then

Call FN2g(M, N, X(), 0, YY())

ElseIf IRev = 2 Then

Call JN2g(M, N, X(), 0, YYp())

End If

Loop While IRev <> 0

Debug.Print "C1, C2 =", X(0), X(1)

Debug.Print "Info =", Info

End Sub

実行結果: C1, C2 = 241.084896112856 5.44942234058364E-04

Info = 0