最小二乘法算法程序
第3章曲线拟合的最小二乘法计算方法
最小二乘拟合,特别是多项式拟合,是最流行的数据处理 方法之一.它常用于把实验数据(离散的数据)归纳总结为经 验公式(连续的函数),以利于进一步的推演分析或应用.
1
结束
§3.2 线性拟合和二次拟合函数
1. 线性拟合
计 已知数据点为 ( xi , yi ), i 1,2,..., n
算 用直线 p( x) a bx作为近似曲线,均方误差为
计
i xi yi xi yi xi2 xi2yi xi3
xi4
0 3 5 15 9 45 27
81
算
1 5 2 10 25 50 125 625
方
2 6 1 6 36 36 216 1296
法
3 8 2 16 64 128 512 4096
课
4 10 4 40 100 400 1000 10000
件
Y ln y, A ln a Y A bx
8
i
xi
0
1
yi
Yi
15.3
2.7279
xi2
xiYi
1
2.7279
1
2
20.5
3.0204
4
6.0408
计
2
3
27.4
3.3105
9
9.9315
算
3
4
36.6
3.6000
16
14.4000
方
4
5
49.1
3.8939
25
19.4695
法
5
6
65.6
4
例1 设5组数据如下表,用一多项式对其进行拟合。
x 3 5 6 8 10
计
最小二乘修正算法
Y Xθ ε
T T T
——————④
于是得：
θ ( X) Y ( X) ε
T T T T
1
1
———⑤ ———⑥
取上式等号右边第一项作为θ的估计量，记为
ˆ (T X) 1 T Y θ IV
式中，Ω称为辅助变量矩阵，Ω中的元素称为辅助变量。
式⑥在形式上与一般最小二乘估计公式相似，但 它在{ε(k)}为有色噪声情况下是参数的一致无偏估计量。 这是因为式⑤有
ˆ
IV
( X)
T 1 T
———⑦
利用滤波法的出发点是：如果辅助模型是未 受噪声ξ(k)干扰的原系统模型，那么辅助模型的 输出就是v(k),且满足： 1 显然，由{v(k)}，{u(k)}构成的下列辅助矩 阵
T (n 1) v(n) v(1) u ( n) u (1) T (n N ) v(n N 1) v( N ) u (n N 1) u ( N )
B(q ) v(k ) u ε无关而与X相关，因此②，③式 必满足。
但被辨对象的理想输出yu(k)得不到，因此只能利 ˆ ˆ 用被辨系统的一个拟合模型 B(q 1 ) / A(q 1 ) 代替 B(q1 )/ A(q1 ) ，并利用这个辅助模型的输出来代替v(k)。 ˆ (q 1 ) 显然，辅助模型的输出应为v(k ) B ˆ u (k ) 1 ˆ A(q ) 而辅助变量矩阵相应为
A(q 1 ) y(k ) B(q 1 )u(k ) (k )
其向量形式为：
Y X
————①
• 式中
xT (n 1) y (n 1) (n 1) , ε , X Y xT ( n N ) y (n N ) (n N )
最小二乘法计算方法
最小二乘法计算方法最小二乘法（Least Squares Method）是一种用于拟合数据和求解最优参数的数学方法。
它被广泛应用于各个领域，如物理学、工程学、经济学等。
本文将介绍最小二乘法的基本原理、应用领域以及计算步骤。
最小二乘法的基本原理是通过最小化数据与拟合函数之间的误差平方和来确定最优参数。
对于一个给定的数据集，我们希望找到一个函数，使得该函数与数据之间的误差最小。
最小二乘法的核心思想是，通过调整函数的参数，使得误差平方和达到最小值。
最小二乘法可以应用于各种函数形式的拟合，包括线性函数、多项式函数、指数函数等。
在实际应用中，我们常常使用线性函数进行拟合，因为线性函数的计算较为简单，且可以用来拟合各种数据。
最小二乘法的应用领域非常广泛。
在物理学中，最小二乘法可以用来拟合实验数据，从而获得物理模型的参数。
在工程学中，最小二乘法可以用来优化控制系统的参数，提高系统的性能。
在经济学中，最小二乘法可以用来分析经济数据，预测经济趋势。
下面我们将介绍最小二乘法的计算步骤。
首先，我们需要确定拟合函数的形式。
对于线性函数拟合，拟合函数的形式可以表示为：y = a + bx，其中a和b是待确定的参数。
然后，我们需要收集实验数据，并将数据表示为坐标系中的点。
接下来，我们需要计算每个数据点到拟合函数的垂直距离，并将这些距离的平方求和，得到误差平方和。
最后，我们使用数学方法（如求导）来确定误差平方和的最小值，并得到最优参数a和b。
最小二乘法的计算步骤可以总结为以下几步：1. 确定拟合函数的形式；2. 收集实验数据，并将数据表示为坐标系中的点；3. 计算每个数据点到拟合函数的垂直距离，并求和得到误差平方和；4. 使用数学方法求解误差平方和的最小值，并得到最优参数。
需要注意的是，最小二乘法并不一定能得到唯一的最优解。
在实际应用中，我们需要综合考虑其他因素，如数据的可靠性、拟合函数的合理性等。
最小二乘法作为一种常用的数据拟合和参数求解方法，具有广泛的应用前景。
递推最小二乘法算法
题目(递推最小二乘法)考虑如下系统：y(k) -1.5y(k -1) 0.7y(k -2) = u(k -3) 0.5u(k -4) (k)式中，(k)为方差为0.1的白噪声。
取初值P(0) =1061、二(0) = 0。
选择方差为1的白噪声作为输入信号u(k)，米用PLS法进行参数估计。
Matlab代码如下：clear allclose allL=400; %仿真长度uk=zeros(4,1); %输入初值：uk(i)表示u(k-i) yk=zeros(2,1); %输出初值u=randn (L,1); %输入采用白噪声序列xi=sqrt(0.1)*randn(L,1); % 方差为0.1 的白噪声序列theta=[-1.5;0.7;1.0;0.5]; %对象参数真值thetae_ 1= zeros(4,1); % ()初值P=10A6*eye(4); %题目要求的初值for k=1:Lphi=[-yk;uk(3:4)]; %400 4矩陳phi第k行对应的y(k-1),y(k-2),u(k-3), u(k-4) y(k)=phi'*theta+xi(k); % 采集输出数据%递推最小二乘法的递推公式K=P*phi/(1+phi'*P*phi);thetae(:,k)=thetae_1+K*(y(k)-phi'*thetae_1);P=(eye(4)-K*phi')*P;%更新数据thetae_1=thetae(:,k);for i=4:-1:2uk(i)=uk(i-1);enduk(1)=u(k);for i=2:-1:2火i)=yk(i-1);end yk(1)=y(k);endplot([1:L],thetae); %line([1,L],[theta,theta]); xlabel('k'); ylabel('参数估计a、b'); lege nd('a_1','a_2','b_0','b_1');%axis([0 L -2 2]);结果分析如下系统方程为y(k) -1.5y(k -1) 0.7y(k -2) = u(k -3) 0.5u(k -4) (k)由CAR模型：y(k) =1.5y(k _1) _0.7y(k _2) u(k _3) 0.5u(k _4) (k)可以得到：-1-1 -2A(z ) =1 1.5z -0.7zB(z 斗)=z，(1 0.5z4)我们可以知道纯迟延为3T，na=2，nb=1，d=3，phi(k,:)=[-yk;uk(3:4)]'; 400 4矩阵phi第k行对应的y(k-1),y(k-2),u(k-3), u(k-4)y(k)=phi*theta+xi(k);输出等于矩阵phi与真值矩阵相乘加上白噪声。
最小二乘法分类
最小二乘法分类最小二乘法（Least Squares Method）是一种常用的参数估计方法，用于寻找一个函数模型的最佳拟合参数，使得模型的预测值与观测值的残差平方和最小化。
这种方法最早由高斯提出，并被广泛应用于统计学和计算机科学等领域。
本文将介绍最小二乘法的基本原理、应用场景以及相关的算法和评估指标。
一、基本原理：最小二乘法用于求解形如y = f(x;θ) 的函数模型的参数θ，其中y是观测值，x是自变量，f是函数模型。
最小二乘法的目标是找到最佳的参数θ，使得模型的预测值与实际观测值之间的残差平方和最小化。
具体步骤如下：1. 定义函数模型：根据具体问题，选择适当的函数模型，如线性模型、多项式模型、指数模型等。
2. 表达目标函数：根据函数模型和参数θ，将目标函数表达为关于θ的函数形式。
3. 定义损失函数：通常采用残差的平方和作为损失函数，即Loss = Σ(y_i - f(x_i;θ))^2 。
4. 求解参数θ：通过最小化损失函数，即求解使得∂Loss/∂θ = 0 的参数θ。
5. 参数估计：根据求解得到的参数θ，即可获得最佳的函数模型。
二、应用场景：最小二乘法在各个领域都有广泛的应用，以下是一些常见的应用场景：1. 线性回归：最小二乘法用于拟合线性回归模型，求解自变量与因变量之间的关系。
2. 特征选择：最小二乘法可用于特征选择，筛选对目标变量影响最大的特征。
3. 数据压缩：通过最小二乘法可以估计出一个低维子空间，将高维数据进行压缩。
4. 图像处理：最小二乘法可用于图像去噪、图像恢复等问题，如使用低秩矩阵模型对图像进行恢复。
5. 信号处理：最小二乘法可用于信号滤波、信号恢复等问题，如基于 DCT 的音频和图像压缩。
三、算法与评估指标：1. 最小二乘法的数值解：在实际应用中，最小二乘法的数值解可以通过各种数值优化算法来求解，包括梯度下降法、牛顿法、共轭梯度法等。
2. 算法评估指标：常用的评估指标包括残差平方和（Residual Sum of Squares, RSS）、均方误差（Mean Square Error, MSE）以及决定系数（Coefficient of Determination, R^2）等。
最小二乘法
y (k ) = x(k ) + v(k )
y (k ) +
x(k ) = y (k ) − v(k )
n i=0
∑
n
i =1
ai y (k − i) =
∑
biu ( k − i ) + v ( k ) +
∑
n
i =1
aiv (k − i)
ξ (k ) = v(k ) +
∑
n
i =1
a iv (k − i)
系统辨识——最小二 (n) L − y (1) u (n + 1) L y (n + 1) y (n + 2) − y (n + 1) L − y (2) u (n + 2) L = M M M M M y (n + N ) − y (n + N − 1) L − y ( N ) u (n + N ) L a1 a 2 u (1) M ξ (n + 1) u (2) an ξ (n + 2) + b0 M u ( N ) b1 ξ (n + N ) M bn
系统辨识——最小二乘算法 系统辨识——最小二乘算法
系统参数a 系统参数a
a 1.5 a1 a2
^
1
a1 = −1.5003
a2 = 0.7006
^
0.5
a1=-1.5, a2= 0.7
0 theta -0.5 -1 -1.5 -2 0
10
20
30 time steps
40
50
几种最小二乘法递推算法的小结
递推最小二乘法的一般步骤：1. 根据输入输出序列列出最小二乘法估计的观测矩阵ϕ：] )(u ... )1( )( ... )1([)(T b q n k k u n k y k y k ------=ϕ没有给出输出序列的还要先算出输出序列。
本例中， 2)]-u(k 1),-u(k 2),-1),-y(k -[-y(k )(T =k ϕ。
2. 给辨识参数θ和协方差阵P 赋初值。
一般取0θ=0或者极小的数，取σσ,20I P =特别大，本例中取σ=100。
3. 按照下式计算增益矩阵G ：)()1()(1)()1()(k k P k k k P k G T ϕϕϕ-+-= 4. 按照下式计算要辨识的参数θ：)]1(ˆ)()()[()1(ˆ)(ˆ--+-=k k k y k G k k T θϕθθ5. 按照下式计算新的协方差阵P ：)1()()()1()(---=k P k k G k P k P T ϕ6. 计算辨识参数的相对变化量，看是否满足停机准则。
如满足，则不再递推；如不满足，则从第三步开始进行下一次地推，直至满足要求为止。
停机准则：εϑϑϑ<--)(ˆ)1(ˆ)(ˆmax k k k i i i i 本例中由于递推次数只有三十次，故不需要停机准则。
7. 分离参数：将a 1….a na b 1….b nb 从辨识参数θ中分离出来。
8. 画出被辨识参数θ的各次递推估计值图形。
为了说明噪声对递推最小二乘法结果的影响，程序5-7-2在计算模拟观测值时不加噪声， 辨识结果为a1 =1.6417，a2 = 0.7148，b1 = 0.3900，b2 =0.3499，与真实值a1 =1.642， a2 = 0.715， b1 = 0.3900，b2 =0.35相差无几。
程序5-7-2-1在计算模拟观测值时加入了均值为0，方差为0.1的白噪声序列，由于噪声的影响，此时的结果为变值，但变化范围较小，现任取一组结果作为辨识结果。
最小二乘法实现源码
三、最小二乘法之多项式拟合算法和源程序：f unction p=naf it(x,y,m)A=zeros(m+1,m+1);f or i=0:mf or j=0:mA(i+1,j+1)=sum(x.^(i+j));endb(i+1)=sum(x.^i.*y);enda=A\b';p=f liplr(a');四、实验用例：假定某天的气温变化记录如下表，试用最小二乘法找出这一天的气温变化规律.t/h 0 1 2 3 4 5 6 7 8 9 10 11 1215 14 14 14 14 15 16 18 20 22 23 25 28T/ Ct/h 13 14 15 16 17 18 19 20 21 22 23 2431 32 31 29 27 25 24 22 20 18 17 16T/ C考虑下列类型函数，计算误差平方和，并作图比较结果：（1）二次函数；（2）三次函数；（3）四次函数；（4）函数.五、实验结果：输入：x=0:24;y=[15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];结果：（1）：p=naf it(x,y,2)p =-0.09362.59438.4157xi=0:0.1:24;y i=polyv al(p,xi);subplot(2,2,1);plot(x,y,'o',xi,y i,'k');title('polyfit');（2）：p=naf it(x,y,3)p =-0.00800.1931-0.102213.2513y i=polyv al(p,xi);subplot(2,2,2);plot(x,y,'o',xi,y i,'k');title('poly interp');（3）：p=naf it(x,y,4)p =0.0009-0.05210.8658-3.525716.6041 y i=polyv al(p,xi);subplot(2,2,3);plot(x,y,'o',xi,y i,'k');title('linear');（4）：p=naf it(x,log(y),2)p =-0.00450.12532.3866y i=exp(polyv al(p,xi));subplot(2,2,4); plot(x,y,'o',xi,y i,'k');title('linear');title('spline');计算误差平方和：（1）：p=naf it(x,y,2)p =-0.09362.59438.4157y1=polyv al(p,x);e1=sum((y-y1).^2)e1=241.2443（2）：p=naf it(x,y,3)p =-0.00800.1931-0.102213.2513y2=polyv al(p,x);e2=sum((y-y1).^2)e2 =241.2443（3）：p=naf it(x,y,4)p =0.0009-0.05210.8658-3.525716.6041 y3=polyv al(p,x);e3=sum((y-y3).^2)e3=36.2837（4）：p=naf it(x,log(y),2)p =-0.00450.12532.3866y1=exp(polyval(p,x));e4=sum((y-y1).^2)e4=178.6060。