« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »

2018年2月 8日 (木)

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
 
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用

« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »

パソコン・インターネット」カテゴリの記事

学問・資格」カテゴリの記事

日記・コラム・つぶやき」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

トラックバック

« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »

最近のコメント

2018年11月
        1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30  
フォト
無料ブログはココログ