14. 非線形最小二乗法

注 – 例題ワークシートをダウンロードすることができます. 1. 例題ワークシートの入手と使用方法 をご覧ください.

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

例題

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

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

  
  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 Ndata = 14, Nparam = 2
Dim Xdata(Ndata - 1) As Double, Ydata(Ndata - 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(Nparam - 1) As Double, Fvec(Ndata - 1) As Double, Tol As Double
    Dim Info As Long, I As Long
    '--- Input data
    For I = 0 To Ndata - 1
        Xdata(I) = Cells(9 + I, 1)
        Ydata(I) = Cells(9 + I, 2)
    Next
    '--- Start calculation
    Tol = 0.0000000001  '** 1e-10
    For I = 0 To 1
        '--- Initial values
        X(0) = Cells(6 + I, 1)
        X(1) = Cells(6 + I, 2)
        '--- Nonlinear least squares
        Call Lmdif1(AddressOf F, Ndata, Nparam, X(), Fvec(), Tol, Info)
        '--- Output result
        Cells(6 + I, 3) = X(0)
        Cells(6 + I, 4) = X(1)
        Cells(6 + I, 5) = Info
    Next
End Sub

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

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

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

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

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

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

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

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

ソルバーを使用した解き方

XLPackソルバーアドインの「非線形最小二乗法」を使って解くこともできます. C列で B8, C8 のパラメータ値を使ってモデル関数値を計算します. 例えば, 式 =B$8*(1-EXP(-C$8*A11)) がC11セルに入力されています. D列では残差を求めています. 例えば, 式 =B11-C11 がD11セルに入力されています. ソルバーはこの残差から最小二乗解を求めます.

ソルバーについては こちら も参照ください.

公開日: 2015/10/11 (更新日: 2022/11/03)