前言:
這一篇是上一篇 FORTH 語言專題的延續, 所以前言幾乎一樣! (哈..)
是因為這樣的, 在工作上, 台灣過去20年來培養了大量的製程工程師. 有製程工程師就會有大量的工程實驗, 或是製程實驗. 有製程實驗就免不了產生大量的實驗數據要分析, 要歸納. 於是, 大家都一定不陌生的是最小平方法的曲線配適. 然後, 每個工程師們都不陌生就是如何用 Excel 來做這個最小平方法的多項式配適.
這邊, 不甘於只是用 Excel 來做這件事啦! 想說自己來寫一個程式來做這個最小平方法的多項式配適, 這樣一定很方便. 如果你熟悉數值分析技術的話會知道, 多階的最小平方法多項式曲線配適, 最後會轉化變成要解一個n階的線性聯立方程組問題. 所以前一篇, 我們已經把求解任意階的聯立方程組的數值程式給搞定了! 這一篇, 就要直接進入核心, 要開始寫這個最小平方法多項式曲線配適的數值程式囉!
一樣, 為了推廣我所非常喜愛的很稀有的 FORTH 程式語言, 這裡的程式一概都用 FORTH 語言來撰寫, 用來展現一下 FORTH 語言優秀的高階性, 跟交談性, 以及比較少看到的, 用 FORTH 來執行浮點運算, 數值分析程式範例. 示範一下 ANSI FORTH 也可以很適合來做各種物理化學, 數學, 統計, 實驗數據的數值分析應用或模擬計算的!
多項式函數:
這個大家高中就學過了, 應該非常不陌生的! 特別是一元二次多項式啦! 連它等於零的解, 根的公式, 大家應該都是倒背如流的!
y=f(x)= b + a*x 一次多項式
y=f(x)= a0 + a1*x + a2*x^2 二次多項式
y=f(x)= a0 + a1*x + a2*x^2 + a3*x^3 三次多項式
...
廣義的呈現式, n 階一元多項式函數. a0, a1, a2 ... ai 這些是每一項的係數. a0 又被稱之為常數項.
這個多項式函數的特色是, 對於一些實驗數據的曲線, 我們可以用一個n階的多項式函數來逼近他們. 於是這些多項式函數的係數, 變成可以代表這些也許是從數千數萬點實驗數據所得出來的, 他們特性的濃縮的數值. 只要利用少少的多項式函數的幾個係數, 就可以代表及預測數千或數十萬實驗數據的特性, 或它們物裡化學反應的特性, 這就是這個一元多項式函數極為重要的原因.
最小平方法:
假如現在有兩個數據點: { (x1, y1), (x2,y2) }, 我們如何衡量它們之間的距離呢?
方法一: 兩點距離公式 distance = sqrt( (x2-x1)^2 + (y2-y1)^2 ), 其數值直接就是兩點間的絕對距離.
方法二: 兩點相減取絕對值 E = | x2 - x1 | + | y2 - y1 |, 其數值可以評估兩點距離的遠近.
方法三: 兩點距離公式的變形 E = (x2-x1)^2 + (y2-y1)^2, 其數值可以評估兩點距離的遠近.
方法一看似最好, 因為直接是兩點間絕對距離的公式, 物理的因次單位還是距離. 但是因為有開根號, 所以計算複雜. 如果只是單純評估遠近的程度就不適合囉.
方法二也不錯, 因為計算簡單, 物理的因次單位還是距離. 但是缺點是絕對值這個函數無法用微積分運算求極值,
方法三介於之間, 計算較方法二複雜, 比方法一簡單, 且優點是可以透過微積分運算求極值, 這個方法也被稱之為平方法.
所以, 假如現在我們有幾萬個數據點 { (x1, y1), (x2, y2), (x3,y3) ... (xN, yN) } N = 幾萬個, 用方法三, 我們來用平方法來評估這幾萬點的數據點, 它們跟一個已知函數f(x)的距離遠近的數值, 可以用下面這個式子來計算!
E = (y1-f(x1))^2 + (y2-f(x2))^2 + (y3-f(x3))^2 + ... + (yN-f(xN))^2)
每個數據點和函數間的距離:
最小平方法多項式曲線配適:
再回到原來的命題, 假如現在我們有幾萬個數據點 { (x1, y1), (x2, y2), (x3,y3) ... (xN, yN)...} N= 幾萬個. 我們現在想找出一個多項式函數 y = f(x) = y=f(x)= a0 + a1*x + a2*x^2 + a3*x^3 + ..., 讓這個多項式函數, 對這幾萬個數據點, 每一點的 "距離" 最小!
每一點的 "距離" 最小?
看圖說故事, 所以這是一個求極小值的問題! 就是剛剛那個, 用平方法來評估這幾萬點的數據點跟多項式函數 y = f(x) 之間的距離, 求極小值! (所以叫 "最小平方法多項式曲線配適", 名稱的由來!)
所以, 就是在求下面這個距離公式的極小值.
其中多項式函數
為了簡化計算, 我們用4階的多項式函數來推導整個最小平方法的公式吧!
4階的多項式函數:
所以對 N 個 (x, y) 的數據點集合
這些N個數據點對4階的多項式函數平方法距離公式 E(a0,a1,a2,a3,a4) 為
這裡所有的數據點 (Xi, Yi) 都是確定的數值, 所以帶入這個 E(a0,a1,a2,a4) 的公式, 假如 (Xi,Yi) 的數值有幾萬個的話, 公式會有幾萬項. 但是因為都是確定的數值. 真正的未知數只有 a0, a1, a2, a3, a4 這五個4階多項式函數的係數待我們去決定.
我們希望能找到一組 (a0, a1, a2, a3, a4) 的係數數值, 來讓 E(a0,a1,a2,a3,a4) 這個平方法的距離最小. 也就是說, 對 E(a0, a1, a2, a3, a4) 求極小值, 這時候這個 (a0, a1, a2, a3, a4) 所代表的多項式函數, 和所有的數據點 (Xi, Yi) 距離最小!
要對 E(a0, a1, a2, a3, a4) 的函數求極值, 微積分告訴我們, 它們的極值發生在一階導數為零的位置! 所以只要 E(a0, a1, a2, a3, a4) 對 a0, a1, a2, a3, a4 分別偏微分後為零, 我們就可以求得5組線性方程式, 滿足這五組線性方程式的 (a0, a1, a2, a3, a4) 的係數, 其 E(a0, a1, a2, a3, a4, a5) 有極(小)值.
所以, E 對 a0 微分求極值
超越函數微分的 chain rule, 外面微分一次, 乘上裡面微分一次
再化簡
把 Summation 放進去, 所以, 第一個極值位置的方程式為
這一條比較特別, 可以再簡化為
同理, E 對 a1 微分求極值, 得到第二條極值方程式
E 對 a2 微分求極值, 得到第三條極值方程式
E 對 a3 微分求極值, 得到第四條極值方程式
E 對 a4 微分求極值, 得到第五條極值方程式
所以完整要解 a0, a1, a2, a3, a4 這個四階多項式的未知係數, 讓其和萬點已知數據點 (Xi, Yi), 最小平方距離有極值(極小值)的線性聯立方程組如下.
所以, 把這個線性方程組解出來. 最小平方法多項式曲線的4階多項式就被決定出來囉!
這裡, 為了怕嚇到讀者, 所以用四階多項式來舉例. 聰明如你, 馬上可以找到規則, 直接看穿這個最小平方法多項式曲線配適的多項式係數聯立方程組矩陣, 如何推廣到 n 階多項式, 來解它的係數 a0, a1, a2... an 囉!
FORTH程式規劃
前面一篇, 我們已經把可以解任意n階線性聯立方程組的程式碼給寫出來囉 (指令: SolveLinearEQ) 所以這裡很簡單囉, 只要寫個程式, 把數據點 (X1,Y1), (X2,Y2), (X3,Y3) ... (XN,YN) 根據上面的公式, 全部算一算, 填入聯立方程組的矩陣裡面. 最後用指令 SolveLinearEQ 來要求求解 a0, a1, a2, ... an 這樣就解出答案囉!
檢視一下這個矩陣元素, 大概就是把所有的數據點的Xi給多次方後加起來, 或乘上Yi 後加起來. 所以所有的數據點是獨立的, 只用一次. 所以算完之後可以不需要儲存跟保留它們的原值. 要考量, 我們的數據點可能是大數據, 可以是有幾萬個 (Xi, Yi). 所以把這幾萬個值先用電腦的記憶體存起來後再開始做運算, 這是非常不智的! 所以最好的方式就是邊輸入邊計算!
而且要考量到, 這數萬個數據點可能是從文字檔案或CSV 檔案所輸入進來的. 這樣, 最好的方式就是射後不理, 每讀一組 (Xi, Yi) 後就立刻計算並更新陣列裡的所有值, 然後立刻忘掉這組已經讀進來的值, 等待下一組值的進入跟計算.
所以語法設計成這樣, 利用一個輸入指令 Input 來從浮點堆疊裡面輸入 Xi, Yi 兩個資料點數字
4.123e0 7.374e0 Input
Input 會即時的把上面陣列的值跟更新計算出來, 放入 InputBuffer 陣列裡面.
這裡另外一個考量是, 常做多項式曲線配適的人一定有的經驗, 常常不知道要用幾階的多項式來做配適比較適合. 所以常常需要每一階的多項式都要配適看看, 最後看結果好不好才知道要用哪一階多項式比較好.
假如這樣, 每改變階數, 就得改變矩陣. 然後幾萬筆的數據就得重新再輸入一次重新計算. 這樣太沒效率囉!
所以兩個矩陣的使用是需要的, 第一個 InputBuffer 矩陣裡存放最大階數(由常數 MaxDegreeOfPoly來定義)的多項式的所有係數. 所以每用 Input 來輸入一筆資料, 都以最大階數MaxDegreeOfPoly的多項式來預估跟計算, 然後更新存放到 InputBuffer 中.
當使用者要用 PolyFit 來配適跟計算任意階的多項式時, 再從 InputBuffer 的主矩陣裡搬移剛好的參數給 Calculating 矩陣, 再呼叫 SolveLinearEQ 將最後解答解出. 這樣就不用一直重新將所有數據輸入後才能計算囉!
所以程式碼如下
: PolyFit ( degree —)
1+
dup MaxDegreeOfPoly > abort" ..." 超過InputBuffer的維數就abort不做!
dup 2 < abort" ..." 低於1的維數就abort不做!
to degreeOfPoly 存入要計算的多項式階數
Trim&MoveMatrix 將指定多項式階數的資料從 InputBuffer搬到 Calculating矩陣
SolveMatrix 求解 Calculating矩陣
.Result 印出指定階數多項式配適後的結果
;
InputBuffer 矩陣的計算
檢視一下最小平方法極小值的線性聯立方程組矩陣的公式可以發現, 它們是有一定規則的!
係數項斜對角是以 Xi, Xi^2, Xi^3, Xi^4 ... 一直到 Xi^2n (n為多項式的階數) 變化著!
所以
(1) 陣列 Row:0, Col:0 填數據點 (XN, YN) 的總數目 N (例如: 有一萬五千組數據, 就填 N = 15000), 有多少組數據就填多少的數目. 所以每輸入一組數據, 就把這個數字加一.
(2) 係數項斜對角是以 Xi, Xi^2, Xi^3, Xi^4 ... 一直到 Xi^2n (n為多項式的階數) 變化著! 所以可以用兩個 do-loop 來組合出它們的變化, 超過陣列範圍的就捨棄.
(3) 常數項欄位是以 Yi, Xi*Yi, Xi^2*Yi, Xi^3*Yi ... 一直到 Xi^n*Yi (n為多項式的階數) 變化著! 所以簡單用個 do-loop 就搞定囉!
程式碼 for (1), (2)
: UpdateCoefficients ( —)
(2)兩個do-loop實現Xi,Xi^2,Xi^3...,Xi^2n
MaxDegreeOfPoly 2* 0 loop1: i = 0...2n
do i 1+ 0 loop2: j = 0...i+1
do i j i - ( row col) (j, i-j)
確認一下row,col 沒超過陣列範圍,有就捨棄,什麼不做
over MaxDegreeOfPoly < over MaxDegreeOfPoly < and
if ( row col) 計算一下 Xn^j的數值, 加入陣列
InputBufferM j Xn^k f+!
else 2drop then
loop
loop
fnum f@ 0 0 InputBufferM f! (1) 陣列(0,0)填資料總筆數.
;
程式碼 for (3)
: UpdateConstantTerms ( --)
(3) do-loop: i = 0...n , 加 Xn^i*Yn 數值進陣列
MaxDegreeOfPoly 0
do i Xn^k*Yn i MaxDegreeOfPoly InputBufferM f+!
loop
;
所以, 整個資料輸入的過程如下:
第一筆資料:
2.34e0 3.62e0 input
第二筆資料:
4.89e0 8.21e0 input
...
第N筆資料:
54.38e0 7.32e0 input
input 程式碼
: Input ( F: x y — )
Yn f! Xn f! 將資料先放進變數 Xn, Yn方便後面運算!
1 num +! 1e0 fnum f+! 將數據總數加一
UpdateCoefficients 將新輸入的數據計算後加入陣列
UpdateConstantTerms 將新輸入的數據計算後加入陣列
;
Calculating 矩陣的搬移
程式的設計, 如前面所說的, 資料是先不管使用者想要用多少階的多項式來計算, 一開始就以 MaxDegreeOfPoly 所定義的最大階數的多項式來計算的, 所有完整的數據先預先儲存在 InputBuffer 這個陣列. 一直到使用者指定了要計算的階數. 真的要解這個聯立方程組了, 才把對應的數據拷貝到 Calculating 這個計算用途的陣列, 然後呼叫 SolveLinearEQ 來將解答解出.
搬移的程式碼如下
: Trim&MoveMatrix ( —)
\ empty calculating matrix
利用 MatrixDim!重設Calculating 陣列維度至最大
MaxDegreeOfPoly [ MaxDegreeOfPoly 1+ ] literal Calculating MatrixDim!
Calculating EmptyMatrix 以最大維度清空Calculating陣列值
\ re-size matrix
利用 MatrixDim!重設Calculating 至正常大小
degreeOfPoly dup 1+ Calculating MatrixDim!
\ move data from InputBuffer to Calculating, ready for solving related linear eq
從InputBuffer 拷貝係數項至Calculating
degreeOfPoly 0
do degreeOfPoly 0
do i j InputBufferM f@
i j Calculating matrix f! ( Coefficients)
loop
loop
從InputBuffer 拷貝常數項至Calculating
degreeOfPoly 0
do i MaxDegreeOfPoly InputBufferM f@
i degreeOfPoly Calculating Matrix f! ( ConstantTerms)
loop
;
最終結果的驗證
就讓我們用書上找來的, 兩次的實驗數據資料, 用 Excel 來驗證一下, 我們的程式碼有沒有錯誤吧!
使用範例 #1
書上找來的, 實驗數據 AISI C1020 冷滾壓鋼在不同的溫度下,所做的撞擊測試.
左邊是溫度, 右邊是樣本在此溫度下經過撞擊後所吸收的撞擊能量.
這裡, 我們利用最小平方法來對此實驗數據配適4階跟5階的多項式曲線, 一個用 Excel, 另一個使用我們寫的程式. 來驗證我們程式碼的正確性與否.
Excel 作圖後, 結果如下
4階多項式: y = 49.2061305 + 1.46923973 X + 0.0103986946 X^2 -7.04053872E-5 X^3 -7.12085781E-7 X^4
5階多項式: y = 49.2061305 + 1.79429557 X + 0.0103986946 X^2 -2.04786014E-4 X^3 -7.12085781E-7 X^4 + 1.04598974E-8 X^5
來操作我們的程式吧
要正式輸入資料, 使用前請打 EmptyAll, 先把所有資料清掉跟歸零, 2 set-precision 這是個非常方便的 ANSI FORTH 指令, 用來控制符點數在螢幕上顯示的位數. 為讓大家有點感覺, 用 .m 分別把 InputBuffer 跟 Calculating 兩個矩陣裡面的值給印出來.
一開始是空的!
用 input 輸入第一筆 (-100, 4.06) 的資料後, InputBuffer 開始累積 MaxDegreeOfPoly (這裡是15) 階多項式最小平方法配適的極值聯立方程組的係數值.
將所有的實驗值輸入完畢後, 這裡再次 .m 一下 InputBuffer, 印出來給大家參考.
這裡也印一下 Calculating 矩陣, 一開始都是零的!
當你下達 4 PolyFit 指令時, 要求以最小平方法配適 4階多項式. 這時候, 程式會先從 InputBuffer 這個已經有 15階多項式值資料的矩陣, 拷貝 4階多項式值的矩陣值, 放入 Calculating 矩陣.
.m 一下 Calculating 矩陣, 可以清楚地看到這個變化. 最後程式會呼叫 SolveLinearEQ 對 Calculating 求解, 最後求得我們所要的4階多項式.
y = 49.2061305 + 1.46923973 X + 0.0103986946 X^2 -0.0000704053872 X^3 - 0.000000712085781 X^4
這個和 Excel 所算出來的, 完全一致!!!
再打入 5 PolyFit 指令, 來試試 5階多項式吧!
結果是
y = 49.2061305 + 1.79429557 X + 0.0103986946 X^2 -0.000204786014 X^3 -0.000000712085781 X^4 + 0.0000000104598974 X^5 一樣, 這個結果和 Excel 的結果一致!!
哭哭! 發現終端機的截圖, 為了能顯示矩陣的內容, 調太大了! 結果反而 BLOG 縮圖後會看不清楚.
所以底下送上只有資料輸入跟計算結果的新終端機的截圖. 這樣就清楚多囉.
使用範例 #2
我們再找一個書上的範例吧, 書上說, 這是個 ill-condition. 所以假如我們程式碼有問題的話, 結果應該會亂七八糟的!
上面的截圖, 我們的程式碼算出來的結果是
2 polyFit
y = 1.0051 + 0.86418 X + 0.84366 X^2
3 polyFit
y = 0.99991 + 1.0141 X + 0.42526 X^2 + 0.27893 X^3
來用 Excel 驗證一下吧,
2 階多項式
3 階多項式
Bingo!! 即使是 ill condition, 不管是我們的程式碼或是 Excel 兩個算出來結果都是一樣的!!
驗證完畢!! 這個我們所寫出來的程式碼是可信賴, 可以使用的囉!! (大心...
原始程式列表:
\
\ Polynomial least squares curve fitting
\
\ Frank Lin 2018/05/15
\
\
\ assume always with 15 degree of polynomial to perform calculation
\ as f(x) = a0 + a1 x + a2 x^2 + ... a15 x^15
\ in-situ data calculations and update to matrix input-buffer while data inputing
\
\ if you want to extend, change MaxDegreeOfPoly as your need
\
FORTH only FORTH also
\
\ Linear EQ solver math library
\
vocabulary LinearSolver
LinearSolver definitions
: DefMatrix \ define matrix
create ( RowNum ColNum --)
2dup , , * floats allot
;
: ColNum@ ( maddr -- ColNum) \ fetch Column no. of a matrix
@
;
: RowNum@ ( maddr -- RowNum) \ fetch Row no. of a matrix
cell+ @
;
: Matrix ( Row Col maddr -- addr') \ calc the address of the data in matrix
>r r@ @ ( Row Col ColNum | maddr)
rot * ( Col ColNum*Row | maddr)
+ floats [ 2 cells ] literal +
r> +
;
: .m ( maddr --) \ print matrix
cr
dup RowNum@ 0
do dup ColNum@ 0
do dup j i rot Matrix f@ fs. space
loop cr
loop drop
;
: EmptyMatrix ( maddr --) \ empty a matrix, fill zero
dup @ ( maddr ColNum)
swap cell+ dup @ ( ColNum maddr' RowNum)
swap cell+ -rot * floats
erase
;
: DummyMatrix \ pointer to point to a matrix
0 value
;
DummyMatrix MatSelected
: RowScalar* ( f.scalar row --)
MatSelected ColNum@ 0
do dup i MatSelected Matrix dup f@ fover f* f!
loop drop
;
: RowUnitizebyCol ( Row Col --)
2dup MatSelected Matrix f@ ( scalar) fdup f0<>
if 1e0 fswap f/ drop RowScalar*
else 2drop then
;
: Row2Row1Scalar*Subtract>Row2 ( Row2 Row1 f.scalar --)
MatSelected ColNum@ 0
do over i MatSelected Matrix dup f@ ( r2 r1 r2adr f.sc f.r2)
fover ( r2 r1 r2adr f.sc f.r2 f.sc)
over i MatSelected Matrix f@ f* ( r2 r1 r2adr f.sc f.r2 f.r1*sc)
f- f!
loop 2drop fdrop
;
: LowerTrianglize
MatSelected RowNum@ 0
do
MatSelected RowNum@ i
do
i j MatSelected Matrix f@ f0<>
if
i j RowUnitizeByCol
i j =
if j j MatSelected Matrix f@ f0= abort" Pivot is ZERO!!"
else i j 1e0 Row2Row1Scalar*Subtract>Row2 then
then
loop
loop
;
: UpperTrianglize
MatSelected RowNum@ 1
do
0 i
do
i j MatSelected Matrix f@ f0<>
if
i j <>
if i j 2dup MatSelected Matrix f@
Row2Row1Scalar*Subtract>Row2
then
then
-1
+loop
loop
;
: RowSwap ( row2 row1 --)
2dup = if 2drop exit then
MatSelected ColNum@ 0
do over i MatSelected Matrix dup f@ ( r2 r1 r2adr f.r2)
over i MatSelected Matrix dup f@ ( r2 r1 r2adr r1adr f.r2 f.r1)
fswap f! f!
loop 2drop
;
: AdjustPivot ( --)
MatSelected RowNum@ 0
do i i i MatSelected Matrix f@ fabs ( ColMax f.max)
( search for max pivot)
MatSelected RowNum@ i
do i j MatSelected Matrix f@ fabs ( ColMax f.max f.new)
fover fover f<
if fswap drop i then
fdrop
loop
( swap with max pivot)
fdrop i RowSwap
loop
;
: SolveLinearEQ ( m.addr --)
to MatSelected
AdjustPivot
LowerTrianglize
UpperTrianglize
;
LinearSolver also
\
\ Polynomial least squares curve fitting
\
vocabulary PolyFitting
PolyFitting definitions
15 constant MaxDegreeOfPoly \ max degree of polynomial for data store & calculation
0 value degreeOfPoly \ degree of polynomial for calculation
MaxDegreeOfPoly dup 1+ DefMatrix InputBuffer \ main matrix for data store
MaxDegreeOfPoly dup 1+ DefMatrix Calculating \ 2nd matrix for calculation
fvariable Xn \ data inputed, x
fvariable Yn \ data inputed, y
fvariable fnum \ number of data inputed, floating number.
variable num \ number of data inputed, integer.
: f+! ( addr F: x) \ floating point version +!
dup f@ f+ f!
;
: MatrixDim! ( RowNum ColNum maddr —)
\ force to change the row & col number of a matrix
tuck ! cell+ !
;
: InputBufferM InputBuffer Matrix ;
: EmptyAll
0 num ! 0e0 fnum f!
Inputbuffer EmptyMatrix
Calculating EmptyMatrix
;
: Xn^k ( k -- F: Xn^k)
Xn f@ s>d d>f f**
;
: Xn^k*Yn ( k — F: Xn^k*Yn)
Xn^k Yn f@ f*
;
: UpdateCoefficients ( —)
MaxDegreeOfPoly 2* 0
do i 1+ 0
do i j i - ( row col)
over MaxDegreeOfPoly < over MaxDegreeOfPoly < and
if ( row col)
InputBufferM j Xn^k f+!
else 2drop then
loop
loop
fnum f@ 0 0 InputBufferM f!
;
: UpdateConstantTerms ( --)
MaxDegreeOfPoly 0
do i Xn^k*Yn i MaxDegreeOfPoly InputBufferM f+!
loop
;
: Input ( F: x y — )
Yn f! Xn f!
1 num +! 1e0 fnum f+!
UpdateCoefficients
UpdateConstantTerms
;
: Trim&MoveMatrix ( —)
\ empty calculating matrix
MaxDegreeOfPoly [ MaxDegreeOfPoly 1+ ] literal Calculating MatrixDim!
Calculating EmptyMatrix
\ re-size matrix
degreeOfPoly dup 1+ Calculating MatrixDim!
\ move data from InputBuffer to Calculating, ready for solving related linear eq
degreeOfPoly 0
do degreeOfPoly 0
do i j InputBufferM f@
i j Calculating matrix f! ( Coefficients)
loop
loop
degreeOfPoly 0
do i MaxDegreeOfPoly InputBufferM f@
i degreeOfPoly Calculating Matrix f! ( ConstantTerms)
loop
;
: SolveMatrix
Calculating SolveLinearEQ
;
: .Result
cr ." y = " 0 degreeOfPoly Calculating Matrix f@ f.
degreeOfPoly 1
do ." + " i degreeOfPoly Calculating Matrix f@ f. ." x^ " i .
loop cr
;
: PolyFit ( degree —)
1+
dup MaxDegreeOfPoly > abort" Over max degree of poly!"
dup 2 < abort" Below min degree of poly!"
to degreeOfPoly
Trim&MoveMatrix \ move data from inputBuffer to calculating, wait for solving
SolveMatrix \ solving normal linear EQ of least squares curve fitting
.Result \ print result
;
PolyFitting also
\
\ Data for testing
\
vocabulary TestingData
TestingData definitions
\
\ test data #1
\
false
[if]
\
\ fit the experiment curve
\
\
\ Temperature Energy
\ (Degree C) (Joules)
\
9 set-precision
EmptyAll
-100e0 4.06e0 input
-75e0 6.78e0 input
-50e0 9.49e0 input
-25e0 16.27e0 input
0e0 40.67e0 input
25e0 97.62e0 input
50e0 146.43e0 input
75e0 151.85e0 input
100e0 162.70e0 input
4 PolyFit
[then]
\
\ test data #2
\
false
[if]
5 set-precision
emptyAll
0e0 1e0 input
0.25e0 1.284e0 input
0.5e0 1.6487e0 input
0.75e0 2.117e0 input
1e0 2.7183e0 input
2 polyFit
[then]
xxx
留言列表
