squre.m

使用最小二乘法的一次性计算公式来估计y(i) = b0 + b1 * x + b2 * x^2 + b3 * x^3 + ...... bn * x^n中b0, b1, ..., bn的值,并显示你和

% y(i) = b0 + b1 * x + b2 * x^2 + b3 * x^3 + ...... bn * x^n
% 用最小二乘法来估计上式中的b0, b1, ..., bn

x   = ( -4 : .1 : 4 )';    % 横坐标,以列向量表示
Num = length( x );         % 横坐标的点数
A   = [ ones( Num, 1 ), x, x.^2, x.^3, x.^4, x.^5, x.^6, x.^7, x.^8, x.^9, x.^10 ];
% A为如下矩阵:
%    1 x(0) x(0)^2 x(0)^3 x(0)^4 x(0)^5 x(0)^6 x(0)^7 x(0)^8 x(0)^9 x(0)^10
%    1 x(1) x(1)^2 x(1)^3 x(1)^4 x(1)^5 x(1)^6 x(1)^7 x(1)^8 x(1)^9 x(1)^10
%    ......  ......
%    1 x(n) x(n)^2 x(n)^3 x(n)^4 x(n)^5 x(n)^6 x(n)^7 x(n)^8 x(n)^9 x(n)^10

% 生成测试用样本点
y0 = sin( x );       
y  = y0 + 0.5 * randn( Num, 1 );   % 以sin(x)为基础生成一系列随机样本点

% 矩阵A全部参与运算估计[b0, b1, ..., bn]'
% 最后得到的多项式最高次数为10次
b10 = inv( A' * A ) * A' * y;                     
figure( 1 ); clf;  
plot( x, y, '.b' ); hold on;       % 显示样本点
yi = A * b10;                      % 根据估计的系数计算多项式
                                   % y(i) = b0 + b1 * x + ... + bn * x^10
disp( sqrt( mean( ( y0 - yi ).^2 ) ) ); % 计算并显示残差
%E  = y0 - yi;
%V  = ( E' * E ) / 2;
%disp( V );
disp( b10' );                      % 显示估计的多项式系数[b0, b1, ..., bn]
plot( x, yi, 'g' );                % 绘制 x 为10次的图形

% 矩阵前8列参与运算估计[b0, b1, ..., bn]'
% 最后得到的多项式最高次数为7次
b8  = inv( A( :, 1 : 8 )' * A( :, 1 : 8 ) ) * A( :, 1 : 8 )' * y;
figure( 2 ); clf;
plot( x, y, '.b' ); hold on;       % 显示样本点
yi = A( :, 1 : 8 ) * b8;           % 根据估计的系数计算多项式
                                   % y(i) = b0 + b1 * x + ... + bn * x^7
disp( sqrt( mean( ( y0 - yi ).^2 ) ) ); % 计算并显示残差
%E  = y0 - yi;
%V  = ( E' * E ) / 2;
%disp( V );
disp( b8');                        % 显示估计的[b0, b1, ..., bn]
plot( x, yi, 'g' );                % 绘制 x 为7次的图形

% 矩阵前6列参与运算估计[b0, b1, ..., bn]'
% 最后得到的多项式最高次数为5次
b6  = inv( A( :, 1 : 6 )' * A( :, 1 : 6 ) ) * A( :, 1 : 6 )' * y;
figure( 3 ); clf;
plot( x, y, '.b' ); hold on;       % 显示样本点
yi = A( :, 1 : 6 ) * b6;           % 根据估计的系数计算多项式
                                   % y(i) = b0 + b1 * x + ... + bn * x^7
disp( sqrt( mean( ( y0 - yi ).^2 ) ) );
%E  = y0 - yi;
%V  = ( E' * E ) / 2;
%disp( V );
disp( b6');                        % 显示估计的[b0, b1, ..., bn]
plot( x, yi, 'g' );                % 绘制 x 为5次的图形

% 矩阵前4列参与运算估计[b0, b1, ..., bn]'
% 最后得到的多项式最高次数为3次
b4  = inv( A( :, 1 : 4 )' * A( :, 1 : 4 ) ) * A( :, 1 : 4 )' * y;
figure( 4 ); clf;
plot( x, y, '.b' ); hold on;       % 显示样本点
yi = A( :, 1 : 4 ) * b4;           % 根据估计的系数计算多项式
                                   % y(i) = b0 + b1 * x + ... + bn * x^7
disp( sqrt( mean( ( y0 - yi ).^2 ) ) );
%E  = y0 - yi;
%V  = ( E' * E ) / 2;
%disp( V );
disp( b4');                        % 显示估计的[b0, b1, ..., bn]
plot( x, yi, 'g' );                % 绘制 x 为3次的图形

% 矩阵前2列参与运算估计[b0, b1, ..., bn]'
% 最后得到的多项式最高次数为1次
b2  = inv( A( :, 1 : 2 )' * A( :, 1 : 2 ) ) * A( :, 1 : 2 )' * y;
figure( 5 ); clf;
plot( x, y, '.b' ); hold on;       % 显示样本点
yi = A( :, 1 : 2 ) * b2;           % 根据估计的系数计算多项式
                                   % y(i) = b0 + b1 * x + ... + bn * x^7
disp( sqrt( mean( ( y0 - yi ).^2 ) ) );
%E  = y0 - yi;
%V  = ( E' * E ) / 2;
%disp( V );
disp( b2');                        % 显示估计的[b0, b1, ..., bn]
plot( x, yi, 'g' );                % 绘制 x 为1次的图形