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の日) »
「パソコン・インターネット」カテゴリの記事
- Pythonでマイクロストリップパッチアンテナの設計ツールを作る(5) これまでをまとめてクラス化・モジュール化した。パッチサイズ、入力インピーダンス、指向性、効率、利得が計算できる。使い方まとめ。(2024.11.05)
- Pythonでマイクロストリップパッチアンテナの設計ツールを作る(4) アンテナ効率(efficiency)を計算する。放射Q(Qr)、誘電体Q(Qd)、導体Q(Qc)から計算できるが今回は”小型・平面アンテナ(電気情報通信学会)"の式を基にした。(2024.10.31)
- 今夜の選挙結果発表に備えて、出口調査の結果から母比率の95%信頼区間を各候補者についてプロットするものをPythonとSciPyのbinomtest、Matplotlibのerrorbarを使って作る。毎年勝手に予想(カシオの高精度計算サイトの自作式とExcel)していたのがとても楽に。(2024.10.27)
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2024年11月4日 IEEE Antennas and Propagation Magazineでアンテナ技術の進化とアレイのレビュー、EMC EuropeでAIをDe-embeddingに応用、I-PEXがシールド付き細線同軸コネクタ発表、iPhone 17からApple製Wi-Fiチップ?Huaweiの5G-Aᴬ(2024.11.04)
- Pythonでマイクロストリップパッチアンテナの設計ツールを作る(5) これまでをまとめてクラス化・モジュール化した。パッチサイズ、入力インピーダンス、指向性、効率、利得が計算できる。使い方まとめ。(2024.11.05)
- 高周波・RFニュース 2024年10月29日 iFixitがiPad mini 7分解・部品特定、NTTが6Gに向けサブテラヘルツで160Gbps、QorvoのSATCOMウェビナー、STMの新バイオセンサ、GSAのFWAウェビナー、6GとAI、Withwaveの50GHzまで使えるB2Bコネクタ、GapWaveとValeoレーダアンテナ協業(2024.10.29)
「日記・コラム・つぶやき」カテゴリの記事
- 高周波・RFニュース 2024年11月4日 IEEE Antennas and Propagation Magazineでアンテナ技術の進化とアレイのレビュー、EMC EuropeでAIをDe-embeddingに応用、I-PEXがシールド付き細線同軸コネクタ発表、iPhone 17からApple製Wi-Fiチップ?Huaweiの5G-Aᴬ(2024.11.04)
- ライカ監修レンズのほぼデジカメのようなSHARP AQUOS R9 Pro、Xiaomi 14 Ultra、ライカではないがデジカメみたいなOppo Find X8の分解写真(2024.10.30)
- Pythonでマイクロストリップパッチアンテナの設計ツールを作る(5) これまでをまとめてクラス化・モジュール化した。パッチサイズ、入力インピーダンス、指向性、効率、利得が計算できる。使い方まとめ。(2024.11.05)
トラックバック
この記事へのトラックバック一覧です: Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その4:フィッティング(線形):
« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »
コメント