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台!インテルが担当。 »
「パソコン・インターネット」カテゴリの記事
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2024年11月29日 NuvotronicsのPolystrataがQorvoのGaNパワーアンプに採用、WürthがBalunのラインアップ拡張、MinewSemiのモジュールにNordicのnRF54Lが採用、AIチップにPCIe 7.0が必要な理由、SIJの特集はThe Road from 1 Gbps-NRZ to 224 Gbps-PAM4(2024.11.29)
- 高周波・RFニュース 2024年11月27日 Ericsson Mobility Report Nov.2024発行、NordicのnRF54L15がRaytacモジュールに搭載、SIMComの5G RedCapモジュール、北京の清華大学でLevel4の自律ネットワークが稼働、AmpleonのDoherty RF Power Transistor、Rohmのシリコンキャパシタ(2024.11.27)
「日記・コラム・つぶやき」カテゴリの記事
- 高周波・RFニュース 2024年11月29日 NuvotronicsのPolystrataがQorvoのGaNパワーアンプに採用、WürthがBalunのラインアップ拡張、MinewSemiのモジュールにNordicのnRF54Lが採用、AIチップにPCIe 7.0が必要な理由、SIJの特集はThe Road from 1 Gbps-NRZ to 224 Gbps-PAM4(2024.11.29)
- 高周波・RFニュース 2024年11月27日 Ericsson Mobility Report Nov.2024発行、NordicのnRF54L15がRaytacモジュールに搭載、SIMComの5G RedCapモジュール、北京の清華大学でLevel4の自律ネットワークが稼働、AmpleonのDoherty RF Power Transistor、Rohmのシリコンキャパシタ(2024.11.27)
コメント
トラックバック
この記事へのトラックバック一覧です: Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その5:フィッティング(非線形):
« [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
投稿: richarddii | 2018年2月 9日 (金) 09時14分
歴史探偵の気分になれる
ウェブ小説「北円堂の秘密」を知ってますか。
北円堂は古都奈良・興福寺の八角円堂です。
グーグルやスマホで「北円堂の秘密」とネット検索するとヒットし、小一時間で読めます。
その1からラストまで無料です。
順に読めば感動に包まれます。
重複、既読ならご免なさい。
お仕事のリフレッシュや脳トレにご利用下さいませ。
物語と観光地が絡むと興味が倍増します。
投稿: omachi | 2018年2月15日 (木) 17時46分