matlab如何求解多项式

多项式回归

在用回归分析方法作数据拟合时，很多情况下很难写出回归函数的解析表达式，例如股票历史价格的拟合，海岸线拟合，地形曲面拟合等。此时，可借助于多项式回归，根据已知的变量观测数据，构造出一个易于计算的多项式函数来描述变量间的不确定性关系，还可以利用该函数计算非数据节点的变量近似值。

1 多项式回归模型

对于可控变量x和随机变量y的m(m>n)次独立的观测(xi,yi),i=1,2,....,m,若y（因变量，也称为响应变量）和x(自变量)之间的回归模型为：

yi=p1*xi^n+p2*xi^n-1+....+pn*xi+pn+1 +a

a~N(o,aa^2), i=1,2,....m

其中p1,p2,.....,pn+1为未知参数，yi为回归多项式模型

.2 多项式回归的MATLAB实现

（1）polyfit函数的用法

MATLAB中提供了polyfit函数，用来做多项式曲线拟合，求解yi=p1*xi^n+p2*xi^n-1+....+pn*xi+pn+1 +a中的未知参数，polyfit函数的调用格式如下：

<1>p=polyfit(x,y,n)

返回n次（阶）多项式回归方程中系数向量的估计值p，这里的p是一个1x(n+1)的行向量，按降幂排列。输入参数x为自变量观测值向量，y为因变量观测值向量，n为正整数，用来指定多项式的阶数。

<2>[p,s]=polyfit(x,y,n)

还返回一个结构体变量s，可作为polyval函数的输入，用来计算预测值及误差的估计值。s有一个normr字段，字段值为残差的模，其值越小，表示拟合越精确。

<3>[p,s,mu]=polyfit(x,y,n)

首先对自变量x进行标准化变换，x=(x-u)/a ;这里u为x的均值，a为x的标准差，然后对y和标准化变换后的x作多项式回归，返回系数向量的估计值p，结构体变量s,已及mu=[u,a]

(2)polyval函数的用法

MATLAB提供了polyval函数，用来根据多项式系数向量计算多项式的值，调用格式如下：

<1>y=polyval(p,x)

计算n次多项式y=p1*x^n+....+pn*x+pn+1在x处的值，输入参数p=[p1,p2,.....,pn]为系数向量，按降幂排列，x为用户指定的自变量取值向量。输出参数y是与x等长的向量。

<2>[y,delta]=polyval(p,x.s)

根据polyfit函数返回的系数向量p和结构体s计算因变量y的预测值，已经误差的标准差的估计值delta。若误差相互独立，服从同方差的正态分布，则[y-delta,y+delta]可作为预测值的50%置信区间。

<3>y=polyval(p,x,[ ],mu) 或[y,delta]=polyval(p,x,s,mu)

首先对自变量x进行标准化变换，然后进行相应的计算，输入参数mu是polyfit函数的第3个输出参数。

（3）poly2sym函数的用法

MATLAB中提供了poly2sym函数，用来把多项式系数向量转为符号多项式，调用格式如下：

<1> r=poly2sym(p)

根据多项式系数向量p生成多形式的符号表达式r。输入参数p是按降幂排列的多项式系数向量

3 多项式回归示例

现有我国2007年1月到2011年11月的食品价格分类指数数据，如下，根据这些数据研究全国食品零售价格分类指数（上年同月=100）和时间的关系。

（1）数据的散点图

用序号表示时间，记为x，用y表示全国食品零售价格分类指数（上年同月=100）。由于x和y均为一维变量，可以先从x和y的散点图上直观地观察他们之间的关系，然后在做进一步的分析。

[Data,Textdata] = xlsread('食品价格.xls'); %读取数据

x = Data(:,1); %提取Data中的第一列，即观测序号

y = Data(:,3); %提取Data中第3列，即价格指数数据

figure;

plot(x,y,'k.','Markersize',15); %绘制x，y的散点图

xlabel('时间序号');

ylabel('食品零售价格分类指数');

散点图表明x和y的非线性趋势比较明显，可以用多项式曲线进行拟合。

（2）四次多项式拟合

假设y关于x的理论回归方程为：

y=p1*x^4+p2*x^3+p3*x^2+p4*x+p5

其中，p1,p2,p3,p4,p5为未知参数，下面调用polyfit函数求解方程中的未知参数，然后调用poly2sym函数显示多项式的符号表达式

[p4,S4] = polyfit(x,y,4)

r = poly2sym(p4);

r = vpa(r,5) %用vpa函数控制r的精度，保留5位有效数字

p4 =

-0.0001 0.0096 -0.3985 5.5635 94.2769

S4 =

R: [5x5 double]

df: 54

normr: 21.0375

r =

- 0.000074268*x^4 + 0.0096077*x^3 - 0.39845*x^2 + 5.5635*x + 94.277 %与p4向量不完全一致是由于舍入误差的原因

（3）更高次多项式拟合

下面调用polyfit函做更高次(>4)多项式拟合，并把多次拟合的残差的模加以对比

[p5,S5] = polyfit(x,y,5);

S5_normr=S5.normr

[p6,S6] = polyfit(x,y,6);

S6_normr=S6.normr

[p7,S7] = polyfit(x,y,7);

S7_normr=S7.normr

[p8,S8] = polyfit(x,y,8);

S8_normr=S8.normr

[p9,S9] = polyfit(x,y,9);

S9_normr=S9.normr

S5_normr =

21.0359

S6_normr =

16.7662

S7_normr =

12.3067

S8_normr =

11.1946

S9_normr =

10.4050

结果表明，随着多项式次数的提高，残差的模呈下降趋势，单纯从拟合的角度来看，拟合精度会随着多项式次数的提高而提高。

（4）拟合效果图

调用polyval函数计算给定自变量x处的因变量y的预测值，来绘制拟合效果图，从效果图上直观的观察出拟合的准确性

figure;

plot(x,y,'k.','Markersize',15); %x和y的散点图

xlabel('时间序号');

ylabel('食品零售价格分类指数');

hold on;

yd4 = polyval(p4,x);%计算4次多项式拟合的预测值

yd6 = polyval(p6,x);%计算6次多项式

yd8 = polyval(p8,x);% 8次多项式

yd9 = polyval(p9,x);% 9次多项式

plot(x,yd4,'r:+'); %绘制4次多项式拟合曲线

plot(x,yd6,'b--s'); % 6

plot(x,yd8,'g-.d'); % 8

plot(x,yd9,'k-p'); % 9

legend('原始散点','4次多项式拟合','6次多项式拟合','8次多项式拟合','9次多项式拟合')