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

N2p_r() Sub N2p_r ( M As  Long, Md As  Long, N As  Long, X() As  Double, Cov() As  Double, Rd() As  Double, Info As  Long, M1 As  Long, M2 As  Long, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Iout As  Long = 0, Optional Info2 As  Long, Optional NFcall As  Long, Optional NFjcall As  Long, Optional Niter As  Long, Optional S As  Double, Optional NFGcal As  Long, Optional NFcov As  Long, Optional NFjcov As  Long, Optional Rtol As  Double = -1, Optional Atol As  Double = -1, Optional Rdreq As  Long = -1, Optional Covreq As  Long = 0, Optional MaxFcall As  Long = -1, Optional MaxIter As  Long = -1, Optional Dtype As  Long = -1, Optional Dfac As  Double = -1, Optional Dtol As  Double = -1, Optional D0 As  Double = -1, Optional Tuner1 As  Double = -1, Optional Xctol As  Double = -1, Optional Xftol As  Double = -1, Optional Lmax0 As  Double = -1, Optional Lmaxs As  Double = -1, Optional Sctol As  Double = -1, Optional Delta0 As  Double = -1, Optional Dltfdc As  Double = -1  ) 非線形最小二乗法 (適応アルゴリズム) (省メモリ版) (リバースコミュニケーション版) 目的 本ルーチンは, ガウス・ニュートン法, レーベンバーグ・マルカート法などを組み合わせ拡張した適応アルゴリズムにより, M個のN変数非線形関数の二乗和の最小点を求める. min Σfi(x1, x2, ..., xn)^2 (ただし, Σは i = 1 〜 M) IRevに従って関数値およびヤコビ行列をユーザーが与える必要がある. 関数値およびヤコビ行列の計算は, 1回で行わずに複数回に分けて部分ごとに行うことができる. 引数 [in] M データ数. (M > 0) [in] Md IRev = 1 による1回の関数評価で計算できる関数値の最大個数. (0 < Md <= M) [in] N パラメータ数. (0 < N <= M) [in,out] X() 配列 X(LX - 1) (LX >= N) [in] 初期近似解. [out] IRev = 0: 求められた解ベクトル.   IRev = 1, 2: 関数値またはヤコビ行列を求める点.   IRev = 3: 最新の近似解. [out] Cov() 配列 Cov(LCov - 1) (LCov >= N(N + 1)/2) 分散・共分散行列. (下三角部分を列優先の圧縮形式で格納) (Rdreq = 1 または 3で, 正常終了時のみ計算される) [out] Rd() 配列 Rd(LRd - 1) (LRd >= M) 回帰診断ベクトル. (Rd(i)はi番目のデータを削除したときの残差二乗和への影響の程度を表す) (Rdreq = 2 または 3で, 正常終了時のみ計算される) [out] Info = 0: 正常終了. (サブコードをInfo2に返す) = -1: パラメータ M の誤り. (M < N) = -2: パラメータ Md の誤り. (Md < 1 または Md > M) = -3: パラメータ N の誤り. (N < 1) = -4: パラメータ X() の誤り. = -5: パラメータ Cov() の誤り. = -6: パラメータ Rd() の誤り. = -10: パラメータ YY() の誤り. = -11: パラメータ YYp() の誤り. = 7: 特異収束. (近傍のヘッセ行列が特異になった) = 8: 誤収束. (誤った点での収束と思われる. 目標精度が小さすぎる可能性がある) = 9: 関数評価回数の最大値を超えた. = 10: 反復回数の最大値を超えた. = 63: X()の初期点においてF(X)を求めることができない. = 65: X()において微分係数を求めることができない. [out] M1 IRev = 1 または 2 で戻ったとき, 設定すべき最初の関数値あるいはヤコビ行列の最初の行の番号. (1 <= M1 <= M) [out] M2 IRev = 1 または 2 で戻ったとき, 設定すべき最後の関数値あるいはヤコビ行列の最後の行の番号. (N2 = min(M, M1+Md-1)) [in] YY() 配列 YYp(LYY - 1) (LYY >= Nd) IRev = 1の場合, 再呼び出し時にX()におけるM1+i-1番目の関数値をYY(i-1)に与えること (i = 1〜M2-N1+1). [in] YYp() 配列 YYp(LYYp1 - 1, LYYp2 - 1) (LYYp1 >= Nd, LYYp2 >= N) IRev = 2の場合, 再呼び出し時にX()におけるヤコビ行列の第M1+i-1行をYYp(i-1,j-1)に与えること (i = 1〜M2-M1+1, j = 1〜N). [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の時には下記処理を行ってから再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: X()におけるM1+i-1番目の関数値を求め YY(i-1)に設定すること (i = 1〜M2-M1+1). YY()以外の変数を変更してはならない. = 2: X()におけるヤコビ行列の第M1+i-1行をを求め YYp(i-1,j-1)に設定すること (i = 1〜M2-M1+1, j = 1〜N). YYp()以外の変数を変更してはならない. = 3: Iout = 1であれば反復ごとにIRev = 3で戻る. 途中結果(X(), NFcall, NFjcall, Niter, Sなど)を出力する. どの変数も変更してはならない. [in] Iout (省略可) 計算の途中結果を出力するかどうか指定する. (省略時 = 0) = 0: 出力しない (IRev = 3では戻らない). = 1: 反復ごとにIRev = 3で戻り, 途中結果を出力する. (上記以外の値であればIout = 0とみなす) [out] Info2 (省略可) Info = 0 のときのサブコード. = 1: X値の収束条件を満たした. = 2: 関数値相対収束条件を満たした. = 3: X値および関数値相対収束条件の両方を満たした. = 4: 関数値絶対収束条件を満たした. [out] NFcall (省略可) 関数評価回数(IRev = 1で戻った回数). (共分散行列計算のための評価を含む) [out] NFjcall (省略可) ヤコビ行列評価回数(IRev = 2で戻った回数). (共分散行列計算のための評価を含む) [out] Niter (省略可) 反復回数. [out] S (省略可) 求められた解ベクトルX()における残差二乗和. [out] NFGcal (省略可) N2gにおいてサブルーチンFとFjに渡されるNf(呼び出しカウンタ). [out] NFcov (省略可) 共分散行列計算のための関数評価回数(IRev = 1で戻った回数). [out] NFjcov (省略可) 共分散行列計算のためのヤコビ行列評価回数(IRev = 2で戻った回数). [in] Rtol (省略可) 関数値の目標相対精度. (Eps <= Rtol <= 0.1) (省略時 = 1e-10) (以下, Epsはマシンイプシロンを表す) (Rtol < Eps または Rtol > 0.1 であれば省略時の既定値とみなす) [in] Atol (省略可) 関数値の目標絶対精度. (省略時 = 1e-20) (Atol < 0 であれば省略時の既定値とみなす) [in] Rdreq (省略可) 共分散行列と回帰診断ベクトルの計算の指定. (省略時 = 3) = 0: 計算しない. = 1: 共分散行列のみ計算. = 2: 回帰診断ベクトルのみ計算. = 3: 共分散行列と回帰診断ベクトルを計算. (それ以外の値であれば省略時の既定値とみなす) [in] Covreq (省略可) 共分散行列の計算方法. (省略時 = 1) = 1, -1: σ * H^(-1) * (J^T * J) * H^(-1) = 2, -2: σ * H^(-1) = 3, -3: σ * (J^T * J) ただし, σ = S / max(1, M-N) (Sは残差二乗和), H はヘッセ行列, J はヤコビ行列である. 符号は計算に使用するユーザー提供の値を示す. < 0: 関数値(IRev = 1)のみを使って計算する. > 0: 関数値(IRev = 1)およびヤコビ行列(IRev = 2)を使って計算する(ヤコビ行列が使えなければ関数値のみを使用する). (それ以外の値であれば省略時の既定値とみなす) [in] MaxFcall (省略可) 関数Fの呼び出し回数の最大値. (省略時 = 200) (MaxFcall <= 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 反復回数の最大値. (省略時 = 150) (MaxIter <= 0 であれば省略時の既定値とみなす) [in] Dtype (省略可) 自動スケーリングの設定. (Dtype = 0, 1 または 2) (省略時 = 1) = 0: 自動スケーリングを行わない. (スケール係数 = 1) = 1: 反復ごとに自動スケーリングを行う. = 2: 1回目の反復のみ自動スケーリングを行い, その後スケール係数を変更しない. (上記以外の値であれば省略時の既定値とみなす) [in] Dfac (省略可) 自動スケーリングの係数. (0 <= Dfac <= 1) (省略時 = 0.6) 自動スケーリングではD(i)をスケール係数として, すべてのiについてD(i)*X(i)が同程度の大きさになるように反復ごとにD(i)を調整する. まず, D1(i) = max(||Ji||, Dfac*D(i))とする(||Ji||はヤコビ行列行列のi列の2-ノルム). 次に, 以下のようにD(i)を調整する.   D1(i) >= Dtolの場合: D(i) = D1(i)   D1(i) < Dtolの場合: D(i) = D0 (Dfac < 0 または Dfac > 1 であれば省略時の既定値とみなす) [in] Dtol (省略可) 自動スケーリングのしきい値. (Dtol > 0) (省略時 = 1.0e-6) (Dtol <= 0 であれば省略時の既定値とみなす) [in] D0 (省略可) 自動スケーリングの初期値. (D0 > 0) (省略時 = 1) (D0 <= 0 であれば省略時の既定値とみなす) [in] Tuner1 (省略可) 誤収束の判定パラメータ. (0 <= Tuner1 <= 0.5) (省略時 = 0.1) (Tuner1 < 0 または Tuner1 > 0.5 であれば省略時の既定値とみなす) [in] Xctol (省略可) X値の収束判定しきい値. (0 <= Xctol <= 1) (省略時 = Eps^(1/2)) (Xctol < 0 または Xctol > 1 であれば省略時の既定値とみなす) [in] Xftol (省略可) 誤収束の判定しきい値. (0 <= Xftol <= 1) (省略時 = 100*Eps) (Xftol < 0 または Xftol > 1 であれば省略時の既定値とみなす) [in] Lmax0 (省略可) スケーリングされた一番始めのステップ長の最大幅. (Lmax0 > 0) (省略時 = 1) (Lmax0 <= 0 であれば省略時の既定値とみなす) [in] Lmaxs (省略可) [in] Sctol (省略可) LmaxsおよびSctolは特異収束の判定パラメータである. (Lmaxs > 0) (0 <= Sctol <= 1) (省略時: Lmaxs = 1, Sctol = 1e-10) ステップ長 Lmaxs における関数値の減少の推定値が Sctol*abs(f) より小さければ Info = 7 で戻る (fは現在の反復の出発点における関数値). (Lmaxs <= 0 であれば省略時の既定値とみなす) (Sctol < 0 であれば省略時の既定値とみなす) [in] Delta0 (省略可) 共分散行列を計算する際の有限差分ステップサイズを次のようにする (Covreq = 1, 2 の場合). (Eps <= Delta0 <= 1) (省略時 = Sqrt(Eps))   Delta0*max(|X(I)|, 1/D(I))*sign(X(I)) (Delta0 < Eps または Delta0 > 1 であれば省略時の既定値とみなす) [in] Dltfdc (省略可) 共分散行列を計算する際の有限差分ステップサイズを次のようにする (Covreq = -1, -2の場合). (Eps <= Dltfdc <= 1) (省略時 = Eps^(1/3))   Dltfdc*max(|X(I)|, 1/D(I)) (Dltfdc < Eps または Dltfdc > 1 であれば省略時の既定値とみなす) 出典 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 FN2p(M As Long, Md1 As Long, M1 As Long, M2 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 M2 - M1 Fvec(I) = Ydata(M1 + I - 1) - X(0) * (1 - Exp(-Xdata(M1 + I - 1) * X(1))) Next End Sub Sub JN2p(M As Long, Md1 As Long, M1 As Long, M2 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 M2 - M1 Fjac(I, 0) = Exp(-Xdata(M1 + I - 1) * X(1)) - 1 Fjac(I, 1) = -Xdata(M1 + I - 1) * X(0) * Exp(-X(1) * Xdata(M1 + I - 1)) Next End Sub Sub Ex_N2p_r() Const M = 4, Md = 2, N = 2 Dim X(N - 1) As Double, Cov(N * (N + 1) / 2 - 1) As Double, Rd(M - 1) As Double Dim Info As Long, M1 As Long, M2 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 N2p_r(M, Md, N, X(), Cov(), Rd(), Info, M1, M2, YY(), YYp(), IRev) If IRev = 1 Then Call FN2p(M, Md, M1, M2, N, X(), 0, YY()) ElseIf IRev = 2 Then Call JN2p(M, Md, M1, M2, N, X(), 0, YYp()) End If Loop While IRev <> 0 Debug.Print "C1, C2 =", X(0), X(1) Debug.Print "Cov =" Debug.Print Cov(0), Cov(1) Debug.Print Cov(1), Cov(2) Debug.Print "Info =", Info End Sub 実行結果 C1, C2 = 241.084896112857 5.44942234058363E-04 Cov = 20.655501584427 -5.54227857698491E-05 -5.54227857698491E-05 1.49176077119591E-10 Info = 0