统计回归与Matlab软件实现下(非线性回归，多项式回归)

逐步回归

关于变量选择

在有多个自变量的情况下，基于自变量的不同组合可以得到许多回归方程，这些回归方程的效果有好有坏
要得到最优的回归方程

• 回归效果最佳
• 自变量个数尽量少
  变量选择常用方法
• 所有子集回归
• 逐步引入
• 逐步剔除
• 逐步回归

逐步回归的基本思想

基本原则

1. 按照自变量对因变量影响的显著程度，从大到小逐个引入回归方程
2. 每一个变量引入以后，判断先前引入的变量是否由于新变量的引入而变得不显著，若是，则将其剔除
3. 引入一个自变量或剔除一个自变量，为逐步回归的一步
4. 每一步都要进行检验，即，要进行引入变量是否显著以及剔除变量是否不显著的检验分析，一般地：$显著性水平{\alpha }_{进}=0.05，{\alpha }_{出}=0.1$
5. 这个过程反复进行，直至既无变量需要引入，也无变量需要剔除，得到一个最佳的变量组合为止

逐步回归的MATLAB实现

stepwise(x, y, inmodel, alpha)

x，自变量数据，nxm阶矩阵
y，因变量数据，nx1阶矩阵
inmodel，初始模型中包含的自变量子集(缺省时默认为空集)
alpha(缺省时默认为0.05)

逐步回归实例

Matlab程序实现

1. 数据输入

x1=[7 1 11 11 7 11 3 1 2 21 1 11 10]';
x2=[26 29 56 31 52 55 71 31 54 47 40 66 68]';
x3=[6 15 8 8 6 9 17 22 18 4 23 9 8]';
x4=[60 52 20 47 33 22 6 44 22 26 34 12 12]';
y=[78.5 74.3 104.3 87.6 95.9 109.2 102.7 72.5 93.1 115.9 83.8 113.3 109.4]';
x=[x1 x2 x3 x4];

2. 逐步回归

stepwise(x, y)

第一个窗口，表示自变量$x_1,x_2,x_3,x_4$的显著性程度的窗口

• 置信线都是红色的，表示这个变量当前不在模型里面，否则是蓝色，说明当前的模型中没有任何一个自变量，都是常数项

• p-val表示自变量的显著性，p值越小越显著，即$x_4$对因变量的影响最显著

• 看p值是不是小于${\alpha }_{进}=0.05$，小于可以引入
  中间的窗口，表示一些模型的参数和统计量的取值

• 截距，表示常数项
  最后一个窗口，表示逐步回归模型，在历次调试中的RMSE的值，即均方根误差

• 可以直接点击置信线引入

• 或者点击下一步引入
  

再次根据p值判断，引入$x_1$

没有需要剔除的变量
没有需要引入的变量
3. 对变量y与$x_1,x_4$作线性回归

X=[ones(13,1) x1 x4];
[b, bint, r, rint, stats]=regress(y, X)

得

b = 103.0974
    1.4400
    -0.6140

最终模型为
$$y=103.0974+1.4400x_1-0.614x_4$$

可线性化的非线性回归

示例
1.
$$y={\beta }_0+{\beta }_1{e}^{bz}+\epsilon$$
令$x^{\prime }=e^{bx}$
$$y={\beta }_0+{\beta }_1{x}^{\prime }+\epsilon$$
2.
$$y={\beta }_0+{\beta }_{1}x+{\beta }_2{x}^2+\cdots +{\beta }_m{x}^m+\epsilon$$
令$x_1=x,x_2=x^2,\dots ,x_m=x^m$
$$y={\beta }_0+{\beta }_1{x}_1+{\beta }_2{x}_2+\cdots +{\beta }_m{x}^m+\epsilon$$
3.
$$y=a{e}^{bx}{e}^{\epsilon }$$
等式两边取对数
$$\ln y=\ln a+bx+\epsilon$$
能够线性化的非线性回归模型，称为本质线性回归模型

可线性化的非线性回归模型

钢包问题

制作散点图

可以看出是非线性回归或曲线回归问题
可以基于配曲线的思想进行求解

1. 根据散点图确定须配曲线的类型
2. 估计其中的未知参数

常用备选曲线

1. 双曲线$\frac{1}{y}=a+\frac{b}{x}$
   

2. 幂函数曲线$y=a{x}^b$其中$x>0,a>0$
   

3. 指数曲线$y=a{e}^{bx}$，其中$a>0$
   

4. 倒指数曲线$y=a{e}^{b/x}$其中$a>0$
   

5. 对数曲线$y=a+b\ln x,x>0$
   

6. S型曲线$y=\frac{1}{a+b{e}^{-x}}$
   

求解

结合散点图变化趋势，选配倒指数曲线
$$y=a{e}^{b/x}$$

同时取对数，线性化
$$\ln y=\ln a+b\frac{1}{x}$$
$$y^{\prime }=A+b{x}^{\prime }$$
用线性回归方法求解
得
$$\hat{A}=2.4587,\hat{b}=-1.1107$$
由此
$$\hat{a}=e^{\hat{A}}=11.6789$$
最终非线性回归结果为
$$y=11.6789e^{-1.1107/x}$$

当常用的备选曲线都配不上对，如何解决

• 多项式回归
• 神经网络
• 支持向量机

多项式回归

一元多项式回归

基本概念
变量y关于x的m阶多项式回归模型为
$$y=a_1{x}^m+a_2{x}^{m-1}+\cdots +a_m{x}^1+a_{m+1}+\epsilon$$
其中，m是已知的阶数，$a_i(i=1,2,\dots ,m)$是回归参数，$\epsilon \sim N(0,{\sigma }^2)$称为随机误差

一元多项式回归的Matlab实现

1. 参数估计

[p, S]=polyfit(x, y, m)

$x=(x_1,x_2,\dots ,x_n)$为自变量取值向量
$y=(y_1,y_2,\dots ,y_n)$为因变量取值向量
$m$表示多项式回归模型的阶数
$p=(a_1,a_2,\dots ,a_{m+1})$为多项式系数向量
$S$为一个结构数据，用来估计预测误差
2. 点预测和区间预测

y = polyval (p, x)

求回归多项式在x处的预测值y

[y, delta] = polyconf(p, x, S, alpha)

求回归多项式在x处的预测值y及预测值的显著性为$1-alpha$的置信区间：$y\pm delta$，alpha缺省时为0.05

示例

求解

t = 1/30:1/30:14/30;
s=[11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 ... 113.77 129.54 146.48];
[p, S] = polyfit(t, s, 2)

结果

p = 489.2946 65.8896 9.1329
S = struct:
R:[3 x 3 double]
df:11
normr:0.1157

得回归模型
$$\hat{s}=489.2946t^2+64.8896t+9.1329$$
预测及作图

[y_hat, delta] = polyconf(p, t, S)
plot(t, s, 'k+', t, y_hat, 'r')
hold on
plot(t, y_hat-delta, ':b', t, y_hat+delta, ':b')
legend('观测值','预测值','95%置信区间范围')

另解

化为线性回归模型

t = 1/30:1/30:14/30;
s=[11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 ... 113.77 129.54 146.48];
T=[ones(14,1) t'(t.^2)'];
[b, bint, r, rint, stats] = regress(s', T);

多元二项式回归

rstool (x, y, 'model', alpha)

x，nxm矩阵
y，n维列向量
alpha显著性水平(缺省时默认为0.05)
model

示例

x1=[1000 600 1200 500 300 400 1300 1100 1300 300];
x2=[5 7 6 6 8 7 5 4 3 9];
y=[100 75 80 70 50 65 90 100 110 60]';
x=[x1' x2'];
rstool(x, y, 'purequadratic');

输入1000和6

点导出可以得到细节
或输入

beta, rmse

结果

beta = 110.5313
		0.1464
		-26.5709
		-0.0001
		1.8475
rmse = 
		4.5362

故，纯二次多项式回归模型为
$$y=110.5313+0.1464x_1-26.5709x_2-0.0001x_1^2+1.8495x_2^2$$

一般非线性回归

基本概念

非线性回归模型的一般形式
$$y=f(x;\theta )+\epsilon$$
y是因变量
$x=(x_1,x_2,\dots ,x_m)$是自变量向量
$\theta =({\theta }_1,{\theta }_2,\dots ,{\theta }_p$是未知参数向量
$\epsilon$是随机误差项

非线性回归的参数估计

使得残差平方和
$$Q(\theta )=\sum _{i=1}^n(y_i-f(x_i,\theta ){)}^2$$
达到最小参数$\theta$，称为$\theta$的非线性最小二乘估计

非线性回归的MATLAB实现

1. 回归命令

[beta, r, J]=nlinfit(x, y, 'model', beta0)

x为nxm矩阵
y为n维列向量
model表示用m文件定义的非线性函数
beta0表示回归系数初值
beta表示回归系数的估计值
r表示残差
J表示Jacobian矩阵
2. 点预测和置信区间预测

[Y, DELTA] = nlpredci('model', x, beta, r, J)

• 点预测，nlinfit所得的回归函数在x处的预测值Y
• 区间预测，95%置信水平下的置信区间估计：$Y\pm DELTA$

示例

选配倒指数模型
$$y=a{e}^{b/x}$$
求解

1. 编写m文件，实现非线性回归函数$y=a{e}^{b/x}$

function yhat = volum(beta, x)
yhat = beta(1) * exp(beta(2) ./ x);

2. 输入数据

x = 2:16;
y = [6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76];
beta0 = [8 2]';

3. 求回归系数

[beta, r, J] = nlinfit(x', y', 'volum', beta0)
beta

得

beta = 
		11.6036
		-1.0641

4. 预测及作图

[yhat, delta] = nlpredci('volum', x', beta, r, J);
plot(x, y, 'k+', x, yhat, 'r')
hold on
plot(x', yhat-delta, ':b', x', yhat+delta, ':b')