14. 非線形最小二乗法

非線形最小二乗法のXLPackによる解き方の例を示します。

例題

NIST(アメリカ国立標準技術研究所)のホームページ(http://www.itl.nist.gov/div898/strd/)に掲載されているテストデータを使用します。

次のようなモデル関数を使用します。

    f(x) = c1(1 - exp(-c2x))

データは次のとおりです。

y x
10.07 77.6
14.73 114.9
17.94 141.1
23.93 190.8
29.61 239.9
35.18 289.0
40.02 332.8
44.82 378.4
50.76 434.8
55.05 477.3
61.01 536.8
66.40 593.1
75.47 689.1
81.78 760.0

このデータをプロットすると次のようになります。

14-01

テストデータの初期値は2種類指定されており、(500, 0.0001) と (250, 0.0005) です。示されている解は c1=2.3894212918×102, c2=5.5015643181×10-4 です。

VBAプログラムを使用した解き方 (1)

VBAサブルーチンLmdif1を使ったプログラム例を示します。

Const M As Long = 14, N As Long = 2
Dim Xdata(M - 1) As Double, Ydata(M - 1) As Double

Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Dim I As Long
    For I = 0 To M - 1
        Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))
    Next
End Sub

Sub Start()
    Dim X(N - 1) As Double, Fvec(M - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    '--- Input data
    For I = 0 To M - 1
        Xdata(I) = Cells(10 + I, 1)
        Ydata(I) = Cells(10 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(5 + I, 1)
        X(1) = Cells(5 + I, 2)
        '--- Nonlinear least squares
        Call Lmdif1(AddressOf F, M, N, X(), Fvec(), Tol, Info)
        '--- Output result
        Cells(5 + I, 3) = X(0)
        Cells(5 + I, 4) = X(1)
        Cells(5 + I, 5) = Info
    Next
End Sub

目的関数を外部サブルーチンとして用意し、初期値と要求精度を指定してLmdif1を呼び出します。Lmdif1は、初期値として与えられた近似解を出発点に反復計算を行い、その近傍の最小二乗解を求めます。

このプログラムを実行すると次の結果が得られました。なお、得られたパラメータを使った計算値(オレンジの線)とデータ(青い点)をグラフしてあります。

この例ではどちらの初期値の場合も正しい解が求められましたが、実際に使う際には目的の解にできるだけ近い初期値を与える必要があります。そうしないと、問題によっては収束しにくいことがあります。また、局所的な最小点をいくつも持つ場合があり、目的の解以外の局所最小点に収束してしまうことがあります。

VBAプログラムを使用した解き方 (2)

リバースコミュニケーション版のVBAサブルーチンLmdif1_rを使ったプログラム例を示します。

Const M As Long = 14, N As Long = 2
Dim Xdata(M - 1) As Double, Ydata(M - 1) As Double

Sub F(M As Long, N As Long, X() As Double, Fvec() As Double, IFlag As Long)
    Dim I As Long
    For I = 0 To M - 1
        Fvec(I) = Ydata(I) - X(0) * (1 - Exp(-Xdata(I) * X(1)))
    Next
End Sub

Sub Start()
    Dim X(N - 1) As Double, Fvec(M - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    Dim XX(N - 1) As Double, YY(M - 1) As Double, IRev As Long, IFlag As Long
    '--- Input data
    For I = 0 To M - 1
        Xdata(I) = Cells(10 + I, 1)
        Ydata(I) = Cells(10 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(5 + I, 1)
        X(1) = Cells(5 + I, 2)
        '--- Nonlinear least squares
        IRev = 0
        Do
            Call Lmdif1_r(M, N, X(), Fvec(), Tol, Info, XX(), YY(), IRev)
            If IRev = 1 Or IRev = 2 Then
                IFlag = IRev
                Call F(M, N, XX(), YY(), IFlag)
            End If
        Loop While IRev <> 0
        '--- Output result
        Cells(5 + I, 3) = X(0)
        Cells(5 + I, 4) = X(1)
        Cells(5 + I, 5) = Info
    Next
End Sub

目的関数を外部関数として与えるのではなく、IRev = 1のときにXX()の値を使って関数値を計算しYY()に入れて再度Lmdif1_rを呼び出します。リバースコミュニケーション版の詳細についてはこちらを参照してください。

このプログラムを実行すると上と同じ結果が得られます。

Top