在實際工程中,我們常會遇到這種問題:已知一組點的橫縱座標,需要繪製出一條盡可能逼近這些點的曲線(或直線),以進行進一步進行加工或者分析兩個變數之間的相互關係。而獲取這個曲線方程的過程就是曲線擬合。
首先,我們從曲線擬合的最簡單情況——直線擬合來引入問題。如果待擬合點集近似排列在一條直線上時,我們可以設直線 y=ax+by=ax+by=ax+b 為其擬合方程,係數 a=[a,b]a=[a,b]a=[a,b] 為待求解項,已知:
用矩陣形式表達為: y=x0ay=x_ay=x0a,其中:
要求解a,可在方程兩邊同時左乘 x0x_x0 的逆矩陣,如果它是乙個方陣且非奇異的話。
但是,一般情況下 x0x_x0 連方陣都不是,所以我們在此需要用 x0x_x0 構造乙個方陣,即方程兩邊同時左乘 x0x_x0 的轉置矩陣,得到方程: x0ty=x0tx0ax_^y=x_^x_ax0ty=x0tx0a 。
此時,方程的係數矩陣 x0tx0x_^x_x0tx0 為方陣,所以兩邊同時左乘新係數矩陣 x0tx0x_^x_x0tx0 的逆矩陣,便可求得係數向量a ,即:(x0tx0)−1x0ty=a(x_^x_)^x_^y=a(x0tx0)−1x0ty=a 。
方程a=(x0tx0)−1x0tya=(x_^x_)^x_^ya=(x0tx0)−1x0ty 右邊各部分均已知,所以可直接求解得到擬合直線的方程係數向量a。
當樣本點的分布不為直線時,我們可用多項式曲線擬合,即擬合曲線方程為n階多項式
y=∑i=0naixi=anxn+an−1xn−1+...+a1x+a0y=\sum_^a_ix^i=a_nx^n+a_x^+...+a_1x+a_0y=∑i=0naixi=anxn+an−1xn−1+...+a1x+a0 。
用矩陣形式表示為: y=x0ay=x_0ay=x0a ,其中:
待求解項為係數向量a=[an,an−1,...,a2,a1,a0]ta=[a_n,a_,...,a_2,a_1,a_0]^ta=[an,an−1,...,a2,a1,a0]t。
曲線擬合方程y=x0ay=x_0ay=x0a 的求解方法與上面直線的求解方法一樣,也是在方程y=x0ay=x_0ay=x0a 兩邊同左乘x0x_0x0的轉置矩陣得到: x0ty=x0tx0ax_^y=x_^x_ax0ty=x0tx0a,
再同時在新方程兩邊同時左乘x0tx0x_^x_x0tx0 的逆矩陣,得到:(x0tx0)−1x0ty=a(x_^x_)^x_^y=a(x0tx0)−1x0ty=a
上式左邊各部分均已知,所以可直接求解得擬合曲線方程的係數向量a。
%by [email protected]
clear
clcx=[2,4,5,6,6.8,7.5,9,12,13.3,15];
y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
[~,k]=size(x);
for n=1:9
x0=zeros(n+1,k);
for k0=1:k %構造矩陣x0
for n0=1:n+1
x0(n0,k0)=x(k0)^(n+1-n0);
endend
x=x0';
anss=(x'*x)\x'*y';
for i=1:n+1 %answer矩陣儲存每次求得的方程係數,按列儲存
answer(i,n)=anss(i);
endx0=0:0.01:17;
y0=anss(1)*x0.^n ;%根據求得的係數初始化並構造多項式方程
for num=2:1:n+1
y0=y0+anss(num)*x0.^(n+1-num);
endsubplot(3,3,n)
plot(x,y,'*')
hold on
plot(x0,y0)
endsuptitle('不同次數方程曲線擬合結果,從1到9階')
擬合曲線結果:
可以看出看來,當多項式的階數過小是,曲線並不能很好地反映出樣本點的分布情況;但階數過高時,會出現過擬合的情況。
係數矩陣answer:
在matlab中,也有現成的曲線擬合函式polyfit,其也是基於最小二乘原理實現的,具體用法為:ans=polyfit(x,y,n). 其中x,y為待擬合點的座標向量,n為多項式的階數。
下面**是用polyfit函式來做曲線擬合:
clear
clcx=[2,4,5,6,6.8,7.5,9,12,13.3,15];
[~,k]=size(x);
y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
for n=1:9
anss=polyfit(x,y,n); %用polyfit擬合曲線
for i=1:n+1 %answer矩陣儲存每次求得的方程係數,按列儲存
answer(i,n)=anss(i);
endx0=0:0.01:17;
y0=anss(1)*x0.^n ; %根據求得的係數初始化並構造多項式方程
for num=2:1:n+1
y0=y0+anss(num)*x0.^(n+1-num);
endsubplot(3,3,n)
plot(x,y,'*')
hold on
plot(x0,y0)
endsuptitle('不同次數方程曲線擬合結果,從1到9階')
執行結果:
用polyfit擬合的結果與第乙份**執行的結果基本一樣
最小二乘法曲線擬合
設有如下實驗資料x1 2345 6789 1011 1213 141516y 46.4 88.8 9.22 9.59.7 9.86 1010.2 10.32 10.42 10.5 10.55 10.58 10.60 試用最小二乘法多次 1到5次 多項式曲線擬合以上資料。import numpy as...
最小二乘法的曲線擬合
最小二乘法的曲線擬合 bool cdatadistillview leastdoublemultiplication long px,long py,long m,long n,double result,double warp for i 0 i m i z 0 b 0 1 d1 n p 0 c ...
(C )曲線擬合的最小二乘法
using system using system.collections.generic using system.linq using system.text namespace 數值分析實驗報告 region 曲線擬合的最小二乘法 private static void imput conso...