Excel VBAで複素数演算（一次方程式・FFT他, Numerical Recipes移植）、フィッティング（非線形含む）、ルンゲクッタ8次（DOP853)などが使えるライブラリ その5:フィッティング（非線形）

今回の第5回は非線形関数のフィッティングです。

Levenberg-Marquardt Methodを用います。

サブルーチンの引数は前回とよく似ているので同じところは省略ですが、

Sub mrqmin(x() As Double, y() As Double, sig() As Double, ndata As Integer, a() As Double, ia() As Integer, _

      ma As Integer, covar() As Double, alpha() As Double, ByRef chisq As Double, ByRef alamda As Double)

まず初期値を適切に設定し、alamda<0として呼び出します。これで初期化されます。

ここから繰り返しますが、うまく行っていればchi2は小さく、alamdaは10のオーダーで小さくなります。

逆にうまく行ってなければalamdaは10のオーダーで大きくなります。

ここをうまく判定して、最終的に収束したと判定できれば

alamda=0として最後に呼び出し、終了です。

フィッティングするための関数は、VBAが関数を引数に取れないため、名前をきめうちしたサブルーチンとして実装されています。

Sub mrqfuncs(x As Double, a() As Double, ByRef y As Double, dyda() As Double, na As Integer)

となります。パラメータを含む関数形と、そのパラメータごとの微分（ヤコビアン）が必要となります。

例えば、

Sub mrqfuncs(x As Double, a() As Double, ByRef y As Double, dyda() As Double, na As Integer)

    Dim i As Integer

    Dim fac As Double, ex As Double, arg As Double

    y = 0#

    For i = 1 To na - 1 Step 3

        arg = (x - a(i + 1)) / a(i + 2)

        ex = Exp(-arg * arg)

        fac = a(i) * ex * 2# * arg

        y = y + a(i) * ex

        dyda(i) = ex

        dyda(i + 1) = fac / a(i + 2)

        dyda(i + 2) = fac * arg / a(i + 2)

    Next i

End Sub

ではガウシアンでフィッティングします。

例題：

(見てもらうと分かるように単に繰り返しているだけで、収束判定はあえていれていません。ここがポイントなので自らのアルゴリズムを組み込んでください)

    Dim n As Integer, m As Integer

    Dim i As Integer, j As Integer

    Dim x(101) As Double

    Dim y(101) As Double

    Dim sig(101) As Double, dyda(7) As Double

    Dim a(7) As Double, ia(7) As Integer

    Dim covar(7, 7) As Double, chisq As Double, ochisq As Double

    Dim alpha(7, 7) As Double, alamda As Double

    n = 101

    m = 3

    For i = 1 To n

        x(i) = Worksheets("Sheet1").Cells(i + 2, 23)

        y(i) = Worksheets("Sheet1").Cells(i + 2, 24)

        sig(i) = 1#

    Next i

    For i = 1 To m

        a(i) = 3

        ia(i) = 1

    Next i

    alamda = -1

    Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)

    For i = 1 To 10

        Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)

    Next i

    alamda = 0

    Call mrqmin(x, y, sig, n, a, ia(), m, covar, alpha, chisq, alamda)

    Worksheets("Sheet1").Cells(5, 27) = a(1)

    Worksheets("Sheet1").Cells(6, 27) = a(2)

    Worksheets("Sheet1").Cells(7, 27) = a(3)

    For i = 1 To n

        Call mrqfuncs(x(i), a, y(i), dyda, m)

        Worksheets("Sheet1").Cells(i + 2, 25) = y(i)

    Next i

ライブラリ本体：

「ComplexMath.bas」をダウンロード

またルンゲクッタ8次のDOP853ルーチンとそのドライバ。

「DOP853Lib.bas」をダウンロード

「DriverDOP853Lib.bas」をダウンロード

メルセンヌツイスタ用

「SMFTLib.bas」をダウンロード