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の日) »
「パソコン・インターネット」カテゴリの記事
- Wordleを100回やってみて正解率の分布を見た。1回目で当たる確率をpとして、2回目からはk倍ずつ確率が上がっていく(kp,k²p,k³p…)とすると、私の結果はp=0.0089, k=2.4になった。(2022.05.27)
- 高周波回路シミュレータQucsStudioを使ってみる(その5)マルチポート(4ポートだとs4pみたいな)のSパラメータデータをTouchstoneファイルに出力しようとするとExport to SnPでは2ポートになってしまう。Project→Import Dataからポート数を手で入れないと。(2022.05.24)
- LTspiceで高周波解析(2)伝送線路のマイクロストリップライン素子を作ったのでそれをスタブとして使ってスタブBPF(バンドパスフィルタ)のSパラメータを計算してみる。QucsStudioのツールを使ってフィルタ合成したもの。(2022.05.09)
- LTspiceで高周波解析(1)伝送線路のマイクロストリップライン素子(Microstrip LineのWheelerやHammaerstadモデル)をサブサーキットとして作ってみて、かつSパラメータ解析まで.net文をつけてやってみる。線幅、基板厚、比誘電率などをパラメータで計算できる。(2022.05.06)
- ExcelでLAMBDA関数が使えるようになった(17) Sパラメータの標準フォーマットTouchstoneが3ポート以上で並びが変なので、LAMBDA,LET,SEQUENCE,IFS,INDEXなどを組み合わせてフラットにするTouchstoneFlatten(s,n)を作った。ただし今はs3p, s4p(3ポート,4ポート)のみ。(2022.04.22)
「学問・資格」カテゴリの記事
- Wordleを100回やってみて正解率の分布を見た。1回目で当たる確率をpとして、2回目からはk倍ずつ確率が上がっていく(kp,k²p,k³p…)とすると、私の結果はp=0.0089, k=2.4になった。(2022.05.27)
- 高周波回路シミュレータQucsStudioを使ってみる(その5)マルチポート(4ポートだとs4pみたいな)のSパラメータデータをTouchstoneファイルに出力しようとするとExport to SnPでは2ポートになってしまう。Project→Import Dataからポート数を手で入れないと。(2022.05.24)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2022年5月19日 IEEE Microwave Magazineで磁性体を使わない非可逆デバイス、Microwave Journalは5Gスモールセル、THz、Movanoがミリ波血圧&血糖値計発表、三菱の3Dプリンタ衛星アンテナ、など。(2022.05.19)
- シン・ウルトラマンをIMAXで観てきた。面白かった!物理学者の有岡君がラグランジアン,AdS,プランクブレーンとかしゃべってる!Randall–Sundrumも。マグカップは標準理論+重力で、物理監修は予想通り橋本幸士さんでした!ブラックホールの構造もタイムリー!(2022.05.14)
- 個人の正解確率pとして多数決をとったときの正解確率が、pが0.5に近づくとどのくらいの人数で0.9を超えるか計算してみた。p=0.50005だと日本の人口くらいの多数決。(神とさざなみの密室読んで正当性確率が気になったので)(2022.05.10)
「日記・コラム・つぶやき」カテゴリの記事
- 京都府では免許の更新時、優良運転者講習がオンラインでできるようになった!早速やってみたが、PCでやろうとするとマイナンバーカードをNFCリーダで読まないとだめで断念、iPhoneはChromeはだめでSafariじゃないと駄目とかいろいろ罠が…でも最終的にはよかった。(2022.05.28)
- Wordleを100回やってみて正解率の分布を見た。1回目で当たる確率をpとして、2回目からはk倍ずつ確率が上がっていく(kp,k²p,k³p…)とすると、私の結果はp=0.0089, k=2.4になった。(2022.05.27)
- 高周波回路シミュレータQucsStudioを使ってみる(その5)マルチポート(4ポートだとs4pみたいな)のSパラメータデータをTouchstoneファイルに出力しようとするとExport to SnPでは2ポートになってしまう。Project→Import Dataからポート数を手で入れないと。(2022.05.24)
- 新型コロナウイルス、日本の陽性者数&ワクチン接種者数総計をプロット&中国、韓国、アメリカ、ドイツ、フランス、イギリスの陽性者数もプロット(5/22更新)どの国もじわじわと増加は続いている。マスクをしなくなったり油断してるからかな…(2022.05.22)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2022年5月19日 IEEE Microwave Magazineで磁性体を使わない非可逆デバイス、Microwave Journalは5Gスモールセル、THz、Movanoがミリ波血圧&血糖値計発表、三菱の3Dプリンタ衛星アンテナ、など。(2022.05.19)
トラックバック
この記事へのトラックバック一覧です: Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その4:フィッティング(線形):
« Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その3: ガウス・ジョルダンの消去法(複素数) | トップページ | [0,1]の一様乱数を足していって和が1を超えるまでの回数の平均はネイピア数eになる(昨日2/7はeの日) »
コメント