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台!インテルが担当。 »
「パソコン・インターネット」カテゴリの記事
- Xのタイムラインが何かランダムウォークのガウス分布じゃなくて逆正弦法則のように見えてきている…仕方ないのでJulia言語でモンテカルロシミュレーションで試行回数を変えて逆正弦法則とガウス分布を見る。(2025.07.16)
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- Gemini CLIが使えるようになっていたので早速VSCodeのターミナルから使って、JavaScriptで連立一次方程式を計算するコードを書いてもらった。普通にガウスの消去法で計算するhtmlを作ってくれた。(2025.06.27)
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2025年7月17日 IEEE Microwave Magazineで地上監視レーダ解説など、Microwave Journalでミリ波AiPのOTA測定解説など、5G Americasが5G Advancedホワイトペーパー発行、TDKのAI 微細欠陥検出edgeRX Vision、IEEE Trans. on THz の今月号発行(2025.07.17)
- Xのタイムラインが何かランダムウォークのガウス分布じゃなくて逆正弦法則のように見えてきている…仕方ないのでJulia言語でモンテカルロシミュレーションで試行回数を変えて逆正弦法則とガウス分布を見る。(2025.07.16)
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
「日記・コラム・つぶやき」カテゴリの記事
- 高周波・RFニュース 2025年7月17日 IEEE Microwave Magazineで地上監視レーダ解説など、Microwave Journalでミリ波AiPのOTA測定解説など、5G Americasが5G Advancedホワイトペーパー発行、TDKのAI 微細欠陥検出edgeRX Vision、IEEE Trans. on THz の今月号発行(2025.07.17)
- Xのタイムラインが何かランダムウォークのガウス分布じゃなくて逆正弦法則のように見えてきている…仕方ないのでJulia言語でモンテカルロシミュレーションで試行回数を変えて逆正弦法則とガウス分布を見る。(2025.07.16)
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
コメント
トラックバック
この記事へのトラックバック一覧です: 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分