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の日) »
「パソコン・インターネット」カテゴリの記事
- Visual Basic (VB.NET)でC#用の数値計算ライブラリMath.NET Numericsを使う(1)複素行列を定義して一次方程式や逆行列、行列式などを計算する。(2023.03.24)
- KeysightのADSで位相を±180°に限らずに連続にする関数をよく聞かれるがいつも忘れる…unwrap()だ。PythonのNumPyでもあるので(matlabにもある)いい加減に覚えたい。とりあえずPythonでやってみて記憶する。(2023.03.02)
- 50万人が毎年受ける試験で採用、“謎”のプログラミング言語「DNCL」という記事を見た。大学入試センターで2022年1月に出している仕様(日本語プログラムぽい)と、令和7年度試験の問題作成の方向性,試作問題等で出している仕様(Pythonを日本語にしたっぽい)違うのか…(2023.03.01)
- LC共振回路のモンテカルロシミュレーションでL,Cを一様乱数で振って共振周波数を見る、、、と全然一様じゃない。そりゃそうだ。何を勘違いしてたんだろう…(2023.02.27)
「学問・資格」カテゴリの記事
- Visual Basic (VB.NET)でC#用の数値計算ライブラリMath.NET Numericsを使う(1)複素行列を定義して一次方程式や逆行列、行列式などを計算する。(2023.03.24)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース(3/16) IEEE Microwave Magazineは女性マイクロ波研究者特集、Microwave Journalで285GHz帯で30GHz帯域のOTA測定!QorvoがUWB室内ナビのデモ、STMとonsemiのBluetooth新商品、IDTechExの6Gレポート、など。(2023.03.16)
- KeysightのADSで位相を±180°に限らずに連続にする関数をよく聞かれるがいつも忘れる…unwrap()だ。PythonのNumPyでもあるので(matlabにもある)いい加減に覚えたい。とりあえずPythonでやってみて記憶する。(2023.03.02)
- 50万人が毎年受ける試験で採用、“謎”のプログラミング言語「DNCL」という記事を見た。大学入試センターで2022年1月に出している仕様(日本語プログラムぽい)と、令和7年度試験の問題作成の方向性,試作問題等で出している仕様(Pythonを日本語にしたっぽい)違うのか…(2023.03.01)
「日記・コラム・つぶやき」カテゴリの記事
- Visual Basic (VB.NET)でC#用の数値計算ライブラリMath.NET Numericsを使う(1)複素行列を定義して一次方程式や逆行列、行列式などを計算する。(2023.03.24)
- BingのAIチャットに最大の素数はいくつ?と聞くと、存在しないが今知られている最大の素数は2^89589933-1だと答えた。が、これは引用している記事が誤植していて、実際は2^82589933-1だ...間違った記事を引用すると間違った答えになる。英語で聞くと1つ前のメルセンヌ素数に…(2023.03.17)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース(3/16) IEEE Microwave Magazineは女性マイクロ波研究者特集、Microwave Journalで285GHz帯で30GHz帯域のOTA測定!QorvoがUWB室内ナビのデモ、STMとonsemiのBluetooth新商品、IDTechExの6Gレポート、など。(2023.03.16)
- 今日は円周率の日(3.14)ということでChatGPTとBingのAIチャットに円周率を求める面白い方法を教えてくださいと聞いてみた。ChatGPTが訳の分からない方法をRegenerateするたびに提案してくる...ピザ押し。Bingは検索していることもあってまあこんな感じか。(2023.03.15)
トラックバック
この記事へのトラックバック一覧です: Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その4:フィッティング(線形):
« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »
コメント