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

Lmstr_r() Sub Lmstr_r ( 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, XX() As  Double, YY() As  Double, YYpr() As  Double, IRev 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) IRevに従って関数値およびヤコビ行列をユーザーが与える必要がある. ヤコビ行列は行ごとに計算する. 引数 [in] M 関数の数. (M > 0) [in] N 変数の数. (0 < N <= M) [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 初期近似解 [out] IRev = 0: 求められた最終近似解.   IRev = 3: 最新の近似解.   IRev >= 10: ヤコビ行列を求める点. [out] Fvec() 配列 Fvec(LFvec - 1) (LFvec >= M) IRev = 0: 求められた解ベクトルX()における関数値. IRev = 3: 最新の近似解における関数値. [out] Fjac() 配列 Fjac(LFjac1 - 1, LFjac2 - 1) (LFjac1 >= M, LFjac2 >= N) IRev = 0: 上部の 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に返す) = -1: パラメータ M の誤り. (M < N) = -2: パラメータ N の誤り. (N <= 0) = -3: パラメータ X() の誤り. (配列X()の大きさが不足) = -4: パラメータ Fvec() の誤り. (配列Fvec()の大きさが不足) = -5: パラメータ Fjac() の誤り. (配列Fjac()の大きさが不足) = -6: パラメータ FTol の誤り. (FTol < 0) = -7: パラメータ XTol の誤り. (XTol < 0) = -8: パラメータ GTol の誤り. (GTol < 0) = -9: パラメータ Diag() の誤り. (配列Diag()の大きさが不足, または, Mode = 2 のときにDiag(i) < 0) = -11: パラメータ Ipiv() の誤り. (配列Ipiv()の大きさが不足) = -13: パラメータ XX() の誤り. (配列XX()の大きさが不足) = -14: パラメータ YY() の誤り. (配列YY()の大きさが不足) = -15: パラメータ YYpr() の誤り. (配列YYpr()の大きさが不足) = 1: 関数評価回数(IRev = 1)がMaxfevに達した. = 2: 残差二乗和が減少しなくなった. (FTolが小さすぎる) = 3: 解が改善されなくなった. (XTolが小さすぎる) = 4: Fvecとヤコビ行列の列が計算機イプシロン内で直交した. (GTolが小さすぎる) [out] XX() 配列 XX(LXX - 1) (LXX >= N) IRev = 1の場合, 関数値を求めるべき点を返す. [in] YY() 配列 YY(LYY - 1) (LYY >= M) IRev = 1の場合, 再呼び出し時に関数値fi(XX()) (i = 1〜M)を与えること. [in] YYpr() 配列 YYpr(LYYpr - 1) (LYYpr >= N) IRev >= 10の場合, 再呼び出し時にX()におけるヤコビ行列の第(IRev-10+1)行 (dfi/dxj) (i = IRev-10+1, j = 1〜N)を与えること. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: XX()における関数値を求めYY()に設定すること. YY()以外の変数を変更してはならない. = 3: 途中結果(X(), Fvec()など)を出力する (Nprint > 0 の場合). どの変数も変更してはならない. >= 10: X()におけるヤコビ行列の第(IRev-10+1)行を求め YYpr()に設定すること. YYpr()以外の変数を変更してはならない. [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回の反復ごと, 最終反復後に途中結果の出力のために IRev = 3 として戻る. <= 0: 途中結果の出力のために戻らない. (省略時 = 0) [out] Nfev (省略可) 関数評価回数 (IRev = 1). [out] Njev (省略可) ヤコビ行列の評価回数 (IRev >= 10). 出典 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 FLmstr(M As Long, N As Long, X() As Double, Fvec() As Double, Fjrow() 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 And IFlag <= M + 1 Then Fjrow(0) = Exp(-Xdata(IFlag - 2) * X(1)) - 1 Fjrow(1) = -Xdata(IFlag - 2) * X(0) * Exp(-X(1) * Xdata(IFlag - 2)) End If End Sub Sub Ex_Lmstr_r() 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 Dim XX(N - 1) As Double, YY(M - 1) As Double, YYpr(N - 1) As Double Dim IRev As Long, IFlag As Long FTol = 0.00000001: XTol = 0.00000001: GTol = 0.00000001 '1.0e-8 X(0) = 500: X(1) = 0.0001 Mode = 1 IRev = 0 Do Call Lmstr_r(M, N, X(), Fvec(), Fjac(), FTol, XTol, GTol, Diag(), Mode, Ipvt(), Info, XX(), YY(), YYpr(), IRev) If IRev = 1 Then IFlag = 1 Call FLmstr(M, N, XX(), YY(), YYpr(), IFlag) ElseIf IRev >= 10 Then IFlag = IRev - 10 + 2 Call FLmstr(M, N, XX(), YY(), YYpr(), IFlag) End If Loop While IRev <> 0 Debug.Print "C1, C2 =", X(0), X(1) Debug.Print "Info =", Info End Sub 実行結果 C1, C2 = 241.084896089973 5.44942234119956E-04 Info = 0