« [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) | トップページ | 平昌オリンピック開会式のドローンは1218台!インテルが担当。 »

2018年2月 9日 (金)

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

« [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) | トップページ | 平昌オリンピック開会式のドローンは1218台!インテルが担当。 »

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

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

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

コメント

Hello, my name is Richard Dicks!

I`m an academic writer and I`m going to change your lifes onсe and for all
Writing has been my passion since early years and now I can`t imagine my life without it.
Most of my books were sold throughout Canada, USA, Old England and even Russia. Also I`m working with services that help people to save their time.
People ask me "Sir, Richard Dicks, I need your professional help" and I always accept the request, `cause I know, that only I can save their time!

Professional Writer - Richard Dicks - [url=http://www.tesolteacherfind.com]Tesol Teacher Find[/url] Corps

歴史探偵の気分になれる
ウェブ小説「北円堂の秘密」を知ってますか。
北円堂は古都奈良・興福寺の八角円堂です。
グーグルやスマホで「北円堂の秘密」とネット検索するとヒットし、小一時間で読めます。
その1からラストまで無料です。
順に読めば感動に包まれます。
重複、既読ならご免なさい。
お仕事のリフレッシュや脳トレにご利用下さいませ。
物語と観光地が絡むと興味が倍増します。

コメントを書く

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

トラックバック

« [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) | トップページ | 平昌オリンピック開会式のドローンは1218台!インテルが担当。 »

最近のコメント

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  
フォト
無料ブログはココログ