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

今回の4回目ははフィッティング関数.

(1)　線形　y=a+b*x

(2)　線形    y=a1*f1(x)+a2*f2(x)+・・・+an*fn(x)

(3)  非線形  y=f(a1,a2,･･･, an)

の3種類があります。

まずは線形の(1)はfitというサブルーチンを使います。引数は以下の通り。

Sub fit(x() As Double, y() As Double, ndata As Integer, sig() As Double, mwt As Integer, ByRef a As Double, ByRef b As Double, _

                ByRef siga As Double, ByRef sigb As Double, ByRef chi2 As Double, ByRef q As Double)

データはx(), y()という配列に入っており（全部で ndata個）、

そのデータの標準偏差（yのほう）がsig()に入っているとします。

※ただし、標準偏差がわからないときはmwt=0とすればこのデータを使わずに計算できます。

そのとき、y=a+b*xの形の線形関数にΧ^2を最小化することでフィッティングします。

※パラメータa,bの推定誤差はsiga,sigbに、Χ^2の値はchi2に、当てはまりの良さはqにそれぞれ返ります。

　詳細はNumerical Recipesを参照

事例：

    Dim i As Integer, n As Integer, j As Integer

    Dim x() As Double, y() As Double

    Dim sig() As Double

    Dim a As Double, b As Double

    Dim siga As Double, sigb As Double

    Dim chi2 As Double, q As Double

    n = 7

    ReDim x(n), y(n), sig(n)

    For i = 1 To 7

        x(i) = Worksheets("Sheet1").Cells(i + 1, 6)

        y(i) = Worksheets("Sheet1").Cells(i + 1, 7)

    Next i

    MsgBox x(1)

    Call fit(x, y, n, sig, 0, a, b, siga, sigb, chi2, q)

    Debug.Print a

    Debug.Print b

--

次は一般の線形関数のフィッティングです。

わかりやすいのはこういう多項式フィッティング。

y(x) = a1 + a2*x + a3*x^2 + ・ ・ ・ + aM*x^(M?1)

ですが、より一般的な

y(x)=a1*X1(x)+a2*X2(x)+･･･+aM*XM(x)

の形が可能です。

サブルーチンの引数はこのようになります。

Sub lfit(x() As Double, y() As Double, sig() As Double, ndat As Integer, a() As Double, ia() As Integer, _

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

x() ,y()はそれぞれデータが入った配列、Sigはデータの標準偏差で、

パラメータはa()という配列に入ります(ma個）。y=Σa(i)*fi(x)のような形になります。

ただし、ia()という配列を使って0ならそのパラメータはフィッティングしない、1ならする、のように選ぶことができます。

その他パラメータは共分散行列とカイ二乗の値です。Numerical　Recipesを参照＆あとの事例を参照。

ただし、VBAのサブルーチンが関数を引数に取れない関係上、

フィッティングに使う関数の名前は指定します。ライブラリ中のfuncsというパラメータを

一括で返すサブルーチンを書き換えてください。

多項式の場合は

Sub funcs(x As Double, p() As Double, np As Integer)

    Dim j As Integer

    p(1) = 1#

    For j = 2 To np

        p(j) = p(j - 1) * x

    Next j

End Sub

のようになります。

奇数次のsinで展開するときは、

Sub funcs(x As Double, p() As Double, np As Integer)

    Dim j As Integer

    Dim pi As Double

    pi = 3.14159265358979

    For j = 1 To np

        p(j) = Sin(x * pi * (2 * CDbl(j) - 1))

    Next j

End Sub

とできる。

矩形波を展開する事例：

    Dim i As Integer, j As Integer

    Dim n  As Integer, m As Integer

    Dim x() As Double, y() As Double

    Dim sig() As Double

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

    Dim covar() As Double

    Dim chisq As Double

    Dim p() As Double

    n = 21

    m = 4

    ReDim x(n), y(n), sig(n)

    ReDim a(m), ia(m), p(m)

    ReDim covar(m, m)

    For i = 1 To n

        x(i) = Worksheets("Sheet1").Cells(i + 14, 6)

        y(i) = Worksheets("Sheet1").Cells(i + 14, 7)

        sig(i) = 1#

    Next i

    For i = 1 To m

        a(i) = 0.1

        ia(i) = 1

    Next i

    Call lfit(x, y, sig, n, a, ia, m, covar, chisq)

    For j = 1 To m

        Worksheets("Sheet1").Cells(14 + j, 10) = a(j)

    Next j

    For i = 1 To n

        Call funcs(x(i), p, m)

        y(i) = 0#

        For j = 1 To m

            y(i) = y(i) + p(j) * a(j)

        Next j

        Worksheets("Sheet1").Cells(i + 14, 8) = y(i)

    Next i

ライブラリ本体：

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

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

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

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

メルセンヌツイスタ用

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