matlab中的反三角函数(MATLAB多项式回归)
matlab中的反三角函数(MATLAB多项式回归)(1)polyfit函数的用法.2 多项式回归的MATLAB实现 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为回归多项式模型
多项式回归
在用回归分析方法作数据拟合时,很多情况下很难写出回归函数的解析表达式,例如股票历史价格的拟合,海岸线拟合,地形曲面拟合等。此时,可借助于多项式回归,根据已知的变量观测数据,构造出一个易于计算的多项式函数来描述变量间的不确定性关系,还可以利用该函数计算非数据节点的变量近似值。
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为正整数,用来指定多项式的阶数。
(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次多项式拟合')