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

Lmstr1_r() Sub Lmstr1_r ( M As  Long, N As  Long, X() As  Double, Fvec() As  Double, Fjac() As  Double, Tol As  Double, Ipvt() As  Long, Info As  Long, XX() As  Double, YY() As  Double, YYpr() As  Double, IRev As  Long, Optional Info2 As  Long  ) 非線形最小二乗法 (レーベンバーグ・マルカート法) (省メモリ版) (シンプルドライバ) (リバースコミュニケーション版) 目的 本ルーチンはM個のN変数非線形関数の二乗和の最小点をレーベンバーグ・マルカート法により求める. min Σfi(x1, x2, ..., xn)^2 (ただし, Σは i = 1 〜 M) IRevに従って関数値およびヤコビ行列をユーザーが与える必要がある. ヤコビ行列は行ごとに計算する. Lmstr1_rは, FTol = Tol, XTol = Tol, GTol = 0, Maxfev = 100*(n+1), Mode = 1, Factor = 100, Nprint = 0 としてLmstr_rを呼び出すのに相当する. 引数 [in] M 関数の数. (M > 0) [in] N 変数の数. (0 < N <= M) [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 初期近似解 [out] IRev = 0: 求められた最終近似解.   IRev >= 10: ヤコビ行列を求める点. [out] Fvec() 配列 Fvec(LFvec - 1) (LFvec >= M) IRev = 0: 求められた解ベクトルX()における関数値. [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] Tol 二乗和および解の相対誤差の許容値. (Tol >= 0) [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: パラメータ Tol の誤り. (Tol < 0) = -7: パラメータ Ipiv() の誤り. (配列Ipiv()の大きさが不足) = -9: パラメータ XX() の誤り. (配列XX()の大きさが不足) = -10: パラメータ YY() の誤り. (配列YY()の大きさが不足) = -11: パラメータ 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()以外の変数を変更してはならない. >= 10: X()におけるヤコビ行列の第(IRev-10+1)行を求め YYpr()に設定すること. YYpr()以外の変数を変更してはならない. [out] Info2 (省略可) Info = 0 で戻ったときのサブコード. = 1: 二乗和の相対減少値およびその予測値の両方がFTol以下になった. = 2: 連続する2回の反復の相対誤差がXTol以下になった. = 3: 上の2つ共満たした. 出典 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_Lmstr1_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 Tol As Double, 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 Tol = 0.00000001 '1.0e-8 X(0) = 500: X(1) = 0.0001 IRev = 0 Do Call Lmstr1_r(M, N, X(), Fvec(), Fjac(), Tol, 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