__x?x?niy?;y?__in
(3)求随机误差项方差的估计值.
记ei?yi?yi
为第i个样本观测值的残差.即被解释变量的观测值与估计值之差.则随机误差项方差的估计值为:
^ei????n?2（8）
22证明从略.
至此, 普通最小二乘法一元线性回归模型的参数估计问题得到解决。 M atlab是一种功能强大的系统分析和仿真工 我们选用它作为实现曲线拟合的软件工具 Matlab语言编程实现最小二乘法的思路: (1)输入各参量x、y的测量值(以数组形式 ,这样便于在计算过程中引用);
(2)用M atlab语言中的plot函数x、y的曲线 图,以此图对比典型曲线图,选择合适的经验 ;
(3)按照上例中的方法,选一个系数a,求 b)对它的偏导数,求出其计算表达式;
(4)编写M atlab的M函数,用来完成经验公 待定系数a的计算,该函数输入量为x、y、 量为a、Q,按照由最小二乘法推导出的公式 代入数值由x、y、b计算a、Q;
(5)改变b的取值,多次调用该M函数,比较 结果中的Q值,最小的Q值所对应的a、b值即为 所求。
?改变b的取值#这部分工作也可编一个循
环函数,输入b可能取的区间,计算不同b对应的 Q,再进行比较,保留使Q最小的b及对应的a。 但通常b的改变对Q的影响不是线性的,为方便 观察结果并选择适当的b,?改变b的取值#这部 分工作最好还是编程者自己完成。
最后,只要将得到的函数图像和x、y的曲线 关系图进行对比,就可以直观的看到拟合的效 果了。
另外,Matlab语言提供了一个函数,可以完成 线性曲线拟合,这就是函数polyfi。t函数polyfit 的输入量为x、y、n,其中x、y即为需要建立相互 关系的两个量的测量值,以数组的形式输入,n为 多项式的次数;输出的是多项式系数的行向量,而 得到的多项式是降幂的。
例1
给出一组数据点(xi,yi)列入表1-1中，试用线性最小二乘法求拟合曲线，作出拟合曲线.
表1-1 例1的一组数据(xi,yi)
xi yi -2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6 -192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04 解 （1）在MATLAB工作窗口输入程序
>> x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];
y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04];
plot(x,y,'r*'),
legend('实验数据(xi,yi)') xlabel('x'), ylabel('y'),
title('例7.2.1的数据点(xi,yi)的散点图') 运行后屏幕显示数据的散点图（略）.
（3）编写下列MATLAB程序计算f(x)在(xi,yi)处的函数值，即输入程序
>> syms a1 a2 a3 a4
x=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6]; fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4
运行后屏幕显示关于a1,a2, a3和a4的线性方程组
fi =[ -125/8*a1+25/4*a2-5/2*a3+a4,
-4913/1000*a1+289/100*a2-17/10*a3+a4,
-1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4,
a4, 1/1000*a1+1/100*a2+1/10*a3+a4,
27/8*a1+9/4*a2+3/2*a3+a4, 19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4]
编写构造误差平方和的MATLAB程序
>> y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50
68.04];
fi=[-125/8*a1+25/4*a2-5/2*a3+a4,
-4913/1000*a1+289/100*a2-17/10*a3+a4, -1331/1000*a1+121/100*a2-11/10*a3+a4, -64/125*a1+16/25*a2-4/5*a3+a4, a4, 1/1000*a1+1/100*a2+1/10*a3+a4, 27/8*a1+9/4*a2+3/2*a3+a4,
19683/1000*a1+729/100*a2+27/10*a3+a4, 5832/125*a1+324/25*a2+18/5*a3+a4];
fy=fi-y; fy2=fy.^2; J=sum(fy.^2)
运行后屏幕显示误差平方和如下
J=
(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+289
/100*a2-17/10*a3+a4+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2
为求a1,a2,a3,a4使J达到最小，只需利用极值的必要条件
?J?0 (k?1,2,3,4)，?ak得到关于a1,a2,a3,a4的线性方程组，这可以由下面的MATLAB程序完成，即输入程序
>> syms a1 a2 a3 a4
J=(-125/8*a1+25/4*a2-5/2*a3+a4+1929/10)^2+(-4913/1000*a1+28
9/100*a2-17/10*a3+a4...+171/2)^2+(-1331/1000*a1+121/100*a2-11/10*a3+a4+723/20)^2+(-64/125*a1+16/25*a2-4/5*a3+a4+663/25)^2+(a4+91/10)^2+(1/1000*a1+1/100*a2+1/10*a3+a4+843/100)^2+(27/8*a1+9/4*a2+3/2*a3+a4+328/25)^2+(19683/1000*a1+729/100*a2+27/10*a3+a4-13/2)^2+(5832/125*a1+324/25*a2+18/5*a3+a4-1701/25)^2;
Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3);
Ja4=diff(J,a4);
Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3),
Ja41=simple(Ja4),
运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下
Ja11=
56918107/10000*a1+32097579/25000*a2+1377283/2500*a3+23
667/250*a4-8442429/625
Ja21 =
32097579/25000*a1+1377283/2500*a2+23667/250*a3+67*a4+7
67319/625
Ja31 =
1377283/2500*a1+23667/250*a2+67*a3+18/5*a4-232638/125 Ja41 =
23667/250*a1+67*a2+18/5*a3+18*a4+14859/25
解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序
>>A=[56918107/10000, 32097579/25000, 1377283/2500, 23667/250;
32097579/25000, 1377283/2500, 23667/250, 67; 1377283/2500, 23667/250, 67, 18/5; 23667/250, 67, 18/5, 18];
B=[8442429/625, -767319/625, 232638/125, -14859/25]; C=B/A, f=poly2sym(C)
运行后屏幕显示拟合函数f及其系数C如下
C = 5.0911 -14.1905 6.4102 -8.2574 f=716503695845759/140737488355328*x^3 -7988544102557579/562949953421312*x^2 +1804307491277693/281474976710656*x
-4648521160813215/562949953421312 故所求的拟合曲线为
f(x)?5.0911x3?14.1905x2?6.4102x?8.2574.
（4）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形.输入程序
>> xi=[-2.5 -1.7 -1.1 -0.8 0 0.1 1.5 2.7 3.6];
y=[-192.9 -85.50 -36.15 -26.52 -9.10 -8.43 -13.12 6.50 68.04]; n=length(xi);
f=5.0911.*xi.^3-14.1905.*xi.^2+6.4102.*xi -8.2574; x=-2.5:0.01: 3.6;
F=5.0911.*x.^3-14.1905.*x.^2+6.4102.*x -8.2574;
fy=abs(f-y); fy2=fy.^2; Ew=max(fy), E1=sum(fy)/n, E2=sqrt((sum(fy2))/n)
plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off legend('数据点(xi,yi)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),
title('例7.2.1的数据点(xi,yi)和拟合曲线y=f(x)的图形')
运行后屏幕显示数据(xi,yi)与拟合函数f的最大误差Ew，平均误差E1和均方根误差E2及其数据点(xi,yi)和拟合曲线y=f(x)的图形（略）.
Ew = E1 = E2 =
3.105 4 0.903 4 1.240 9
多重共线性——PLS偏最小二乘法
经典的最小二乘估计，必需满足一些假设条件，多重共线性就是其中的一种。实际上，解释变量间完全不相关的情形是非常少见的，大多数变量都在某种程度上存在着一定的共线性，而存在着共线性会给模型带来许多不确定性的结果。
设回归模型y??0??1x1??2x2????pxp为零的数k0,k1,k2?kp使得k0存在完全共线性,如果k0近似的多重共线性。
??如果矩阵X的列向量存在一组不全
?k1xi1?k2xi2???kpxip?0, i=1,2,…n,则称其
?k1xi1?k2xi2???kpxip?0, i=1,2,…n,则称其存在
PLS(偏最小二乘法)
H.Wold在1975年提出的 偏最小二乘法近年来引起广泛的关注，在解决多重共线性方面能很好的达到目的，偏最小二乘法集中了最小二乘法、主成分分析法和典型相关分析的的优点克服了两种方法的缺点。偏最小二乘法吸取了主成分回归提取主成分的思想，但不同的是主成分回归只是从自变量中去寻找主成分与因变量无关，因而主成分与因变量在算法上关系不密切，从而导致最后主成分在实际应用中无法更好的进一步拟合因变量，偏最小二乘法则是从因变量出发，选择与因变量相关性较强而又能方便运算的自变量的线性组合。 偏最小二乘回归分析 偏最小二乘回归分析 考虑p个因变量L与m个自变量
yp
xm的建模问题。偏最小二
ttxm的线性组合，且尽可能多地
u1并要求u1与t1相关程
乘回归的基本作法是首先在自变量集中提出第一成分1(1是
提取原自变量集中的变异信息）；同时在因变量集中也提取第一成分度达到最大。然后建立因变量
yp与1的回归，如果回归方程已达到满意的精度，则算法中
t止。否则继续第二对成分的提取，直到能达到满意的精度为止。若最终对自变量集提取r个成分r，偏最小二乘回归将通过建立
typ与r的回归式，然后再表示为
typ与原自变量的回
归方程式，即偏最小二乘回归方程式。 偏最小二乘回归MATLAB程序代码
单因变量
function y=pls(pz) [row,col]=size(pz); aver=mean(pz);
stdcov=std(pz); %求均值和标准差 rr=corrcoef(pz); %求相关系数矩阵 úta=zscore(pz); %数据标准化
stdarr = ( pz - aver(ones(row,1),:) )./ stdcov( ones(row,1),:); % 标准化数据结果与zscore()一致
x0=pz(:,1:col-1);y0=pz(:,end); %提取原始的自变量、因变量数据
e0=stdarr(:,1:col-1);f0=stdarr(:,end); %提取标准化后的自变量、因变量数据 num=size(e0,1);%求样本点的个数 temp=eye(col-1);%对角阵 for i=1:col-1
%以下计算 w，w*和 t 的得分向量， w(:,i)= ( e0'* f0 )/ norm( e0'*f0 ); t(:,i)=e0*w(:,i) %计算成分 ti 的得分
alpha(:,i)=e0'*t(:,i)/(t(:,i)'*t(:,i)) %计算 alpha_i ,其中(t(:,i)'*t(:,i))等价于norm(t(:,i))^2
e=e0-t(:,i)*alpha(:,i)' %计算残差矩阵 e0=e; %计算w*矩阵 if i==1
w_star(:,i)=w(:,i); else
for j=1:i-1
最小二乘法在数学模型建立与检验中的运用
