在回归分析中,一般任意的数据都可以用一条曲线来表示,这个曲线可以用某一个高次方的代数多项式 y= a + bx + c(x)2 + …来描述,其中 a , b , …是常数。但是这样通过每个点的曲线是没意义的,也不能表示y和x的真实的相关关系。趋势变化才是两者相关关系的合理解释。在已知y=f(x)的形式时,运用最小二乘法可以计算出它的参数。
现在我们通过求直线回归方程来了解最小二乘法。
假定 y1 为预测值,它和 x之间的关系是 y1 = a + bx,y2 为实测值。根据误差理论,算术平均值是最佳值,由此算出的均方误差最小。我们定义 y3 = y2 - y1为实测值相对于回归直线的误差,那么满足均方误差最小的回归直线方程,就是最佳配置直线。
由上文中的满足均方误差最小的回归直线方程,可以转化为:
∑(y2-bx-a)2 = ∑ (y3)2 = min(表示最小值)
要使等式成立,根据微分学的知识,必须使等式左边对a和b的偏微分分别为零,然后微分得到:
∑(y2-bx-a)x = 0
∑(y2-bx-a) = 0
分解开来:
∑ (xy2) - b∑(x)2 -a∑ x = 0
∑ y2 - b∑ x - Na =0
其中N为数据个数。
解上面的方程组为:
b = (N∑(xy2)-(∑ y2)(∑ x))/(N∑(x)2 - (∑ x)2)
a = ((∑ (x)2)(∑ y2)-(∑ (xy2))(∑ x))/(N∑ (x)2-(∑ x)2)
我们定义 y4 和 x4 为所有 y2 和x的平均值,那么就有
a = y4 - bx4
b = ∑((y2 - y4)*(x - x4))/∑(x - x4)2
举例:采集一段时间静止下的加速度数据,然后求这个拟合曲线。
>> y = load('F:\matlab\impact\s.txt'); %导入加速度数据,这y是一个一列多行的矩阵
>> x = 1:length(y); %获取x的数据序列,x是一个一行多列的矩阵
>> n = length(x); %获取数据长度
>> x = x'; %因为后面要用到矩阵运算,所以x 和y的矩阵维度必须一致,而这个代码就是把x的列转成行
>> x_ave = mean(x); %求x的平均值
>> y_ave = mean(y); %求y的平均值
>> b = sum((y-y_ave).*(x-x_ave))/sum((x-x_ave).^2); %b的计算
>> a = y_ave-b*x_ave; %a的计算
>> r = sum((y-y_ave).*(x-x_ave))/sqrt(sum((x-x_ave).^2).*sum((y-y_ave).^2)); %相关系数计算
>> yhat = a+b*x; %y的估计值
>> lyy = sum((y-y_ave).^2); %总平方和
>> U = sum((yhat-y).^2); %回归平方和
>> Q = lyy-U; %剩余平方和
>> Sy = sqrt(lyy/(n-1)); %总均方和
>> Sq = sqrt(Q); %回归均方差
>> Syhat = sqrt(U/(n-2)); %剩余均方差
>> Sa = Syhat*sqrt(1/n+x_ave^2/sum((x-x_ave).^2)); %a的均方差
>> Sb = Syhat/sqrt(sum((x-x_ave).^2)); %b的均方差
>> xx = min(x):1:max(x); %
>> yy = a+b*xx; %回归直线,拟合直线
>> deltayy = Syhat*sqrt(1+1/n+(xx-x_ave).^2./sum((x-x_ave).^2)); %回归直线误差
>> y1 = yy+deltayy; %误差上限
>> y2 = yy-deltayy; %误差下限
>> plot(xx,yy,'-',x,y,'s');
>> xlabel('X');ylabel('Y');
程序运行得到a = -400.9312 b = 0.0022 那么拟合直线为 y = 0.0022 *x -400.9312
a的值很小,可以认为y的值与x的值没有太大关系。这也符合加速度计静止的事实。