Covar Lmder Lmder1 Lmder1_r Lmder_r Lmdif Lmdif1 Lmdif1_r Lmdif_r Lmstr Lmstr1 Lmstr1_r Lmstr_r N2f N2f1 N2f1_r N2f_r N2g N2g1 N2g1_r N2g_r N2p N2p_r

Lmder() Sub Lmder ( F As  LongPtr, M As  Long, N As  Long, X() As  Double, Fvec() As  Double, Fjac() As  Double, FTol As  Double, XTol As  Double, GTol As  Double, Diag() As  Double, Mode As  Long, Ipvt() As  Long, Info As  Long, Optional Info2 As  Long, Optional Maxfev As  Long = 0, Optional Factor As  Double = 0, Optional Nprint As  Long = 0, Optional Nfev As  Long, Optional Njev As  Long  ) 非線形最小二乗法 (レーベンバーグ・マルカート法) 目的 本ルーチンはM個のN変数非線形関数の二乗和の最小点をレーベンバーグ・マルカート法により求める. min Σfi(x1, x2, ..., xn)^2 (ただし, Σは i = 1 〜 M) 関数値およびヤコビ行列を計算するユーザールーチンが必要である. 引数 [in] F 関数fi(x)の値およびヤコビ行列を求めるユーザーサブルーチンで, 次のように定義すること. Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, Fjac() As Double, IFlag As Long) IFlag = 1の場合: 与えられたX()における関数値fi(X)を求め Fvec(i-1)に設定する(i = 1〜M). それ以外の変数を変更してはならない. IFlag = 2の場合: 与えられたX()におけるヤコビ行列(∂fi/∂xj)を求めFjac()に設定する(i = 1〜M, j = 1〜N). それ以外の変数を変更してはならない. IFlag = 0の場合: Nprintを正の値に設定した場合, 途中結果の出力のためNprint回の反復ごとに呼び出される. 引数に渡された X(), Fvec() などの中間結果を出力する. 渡された変数を変更してはならない. End Sub IFlagの値は実行を強制終了させたい場合以外は変更してはならない. 実行を終了させたい場合には負の値に設定して戻る. [in] M 関数の数. (M > 0) [in] N 変数の数. (0 < N <= M) [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 初期近似解. [out] 求められた解ベクトル. [out] Fvec() 配列 Fvec(LFvec - 1) (LFvec >= M) 求められた解ベクトルX()における関数値. [out] Fjac() 配列 Fjac(LFjac1 - 1, LFjac2 - 1) (LFjac1 >= M, LFjac2 >= N) 上部の N×N 部分に, 対角要素の絶対値が昇順になるように並べ替えを行った上三角行列R が入る.   P^T * (J^T * J)*P = R^T * R ここで, P は並べ替えを表す置換行列, J は最終的に求められたヤコビ行列である. [in] FTol 二乗和の相対誤差の許容値. 二乗和の減少の予測値と実際の値の両方がFTolより小さくなったら終了する. (FTol >= 0) [in] XTol 解の相対誤差の許容値. 続く2回の反復間の相対誤差がXTolより小さくなったら終了する. (XTol >= 0) [in] GTol 関数ベクトルとヤコビ行列の列の直交性の許容値. Fvec()とヤコビ行列のすべての列の角度のcosの絶対値がGTolより小さくなったら終了する. (GTol >= 0) [in,out] Diag() 配列 Diag(LD - 1) (LD >= N) [in] Mode = 2の場合にスケーリング因子を入力する. (Diag(i) > 0) [out] Mode = 1の場合に自動的に設定される. [in] Mode = 1: 変数のスケーリングを自動的に行う. = 2: ユーザーがDiag()に設定した値を使って変数のスケーリングを行う. (他の値であれば Mode = 1 とみなす.) [out] Ipvt() 配列 Ipvt(LIpvt - 1) (LIpvt >= N) 置換行列を定義するピボット情報. [out] Info = 0: 正常終了. (サブコードをInfo2に返す) = -2: パラメータ M の誤り. (M < N) = -3: パラメータ N の誤り. (N <= 0) = -4: パラメータ X() の誤り. (配列X()の大きさが不足) = -5: パラメータ Fvec() の誤り. (配列Fvec()の大きさが不足) = -6: パラメータ Fjac() の誤り. (配列Fjac()の大きさが不足) = -7: パラメータ FTol の誤り. (FTol < 0) = -8: パラメータ XTol の誤り. (XTol < 0) = -9: パラメータ GTol の誤り. (GTol < 0) = -10: パラメータ Diag() の誤り. (配列Diag()の大きさが不足, または, Mode = 2のときにDiag(i) < 0) = -12: パラメータ Ipiv() の誤り. (配列Ipiv()の大きさが不足) = 1: 関数呼び出し(IFlag = 1)回数がMaxfevに達した. = 2: 残差二乗和が減少しなくなった. (FTolが小さすぎる) = 3: 解が改善されなくなった. (XTolが小さすぎる) = 4: Fvecとヤコビ行列の列が計算機イプシロン内で直交した. (GTolが小さすぎる) = 5: ユーザーによる強制終了. (IFlag < 0でFから戻った) [out] Info2 (省略可) Info = 0 で戻ったときのサブコード. = 1: 二乗和の相対減少値およびその予測値の両方がFTol以下になった. = 2: 連続する2回の反復の相対誤差がXTol以下になった. = 3: 上の2つ共満たした. = 4: Fvecとヤコビ行列のすべての列との角度のコサインの絶対値がGTol以下になった. [in] Maxfev (省略可) 関数FのIFlag = 1または2とした呼び出し回数の上限値. (省略時 = 100*(N + 1)) (Maxfev <= 0 であれば省略時の既定値とみなす.) [in] Factor (省略可) 初期ステップの最大値を決める定数. ||Diag×X||2×Factorが0でない場合これを最大値とする. 0の場合Factor自身の値を最大値とする. 通常, 0.1〜100の範囲で, 一般的には100とすればよい. (省略時 = 100) (Factor <= 0 であれば省略時の既定値とみなす.) [in] Nprint (省略可) > 0: 最初の反復, Nprint回の反復ごと, 最終反復後に途中結果の出力のために IFlag = 0としてFを呼び出す. <= 0: 途中結果の出力のためのFの呼び出しを行わない. (省略時 = 0) [out] Nfev (省略可) 関数評価回数 (IFlag = 1としたFの呼び出し回数). [out] Njev (省略可) ヤコビ行列の評価回数 (IFlag = 2としたFの呼び出し回数). 出典 netlib/minpack 使用例 次のデータをモデル関数 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 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_Lmder() Const M = 4, N = 2 Dim X(N - 1) As Double, Fvec(M - 1) As Double, Fjac(M - 1, N - 1) As Double Dim FTol As Double, XTol As Double, GTol As Double, Diag(N - 1) As Double Dim Mode As Long, Ipvt(N - 1) As Long, Info As Long FTol = 0.00000001: XTol = 0.00000001: GTol = 0.00000001 '1.0e-8 X(0) = 500: X(1) = 0.0001 Mode = 1 Call Lmder(AddressOf FLmder, M, N, X(), Fvec(), Fjac(), FTol, XTol, GTol, Diag(), Mode, Ipvt(), Info) Debug.Print "C1, C2 =", X(0), X(1) Debug.Print "Info =", Info End Sub 実行結果 C1, C2 = 241.084896089973 5.44942234119956E-04 Info = 0