数学建模之回归分析加例题详解(MATLAB实现)

一元线性回归

变量之间的关系大致可分为两大类:

  1. 确定性的关系:可以用精确的函数关系来表达。例如矩形面积S与边长a,b的关系。
  2. 非确定性的关系:变量之间既互相联系但又不是完全确定的关系,称为相关关系。例如人的身高与体重、农作物产量与降雨量等的关系。

从数量的角度去研究这种非确定性的关系,是数理统计的一个任务. 包括通过观察和试验数据去判断变量之间有无关系,对其关系大小作数量上的估计、推断和预测,等等.
回归分析就是研究相关关系的一种重要的数理统计方法.

一元正态线性回归模型

只有两个变量的回归分析, 称为一元回归分析;
超过两个变量时称为多元回归分析

变量之间成线性关系时, 称为线性回归;
变量间不具有线性关系时, 称为非线性回归.

若y和x之间大体上呈现线性关系, 可假定 y = a + b x + ϵ y=a+bx+\epsilon y=a+bx+ϵ
其中a 和 b是未知常数, ε表示其它随机因素的影响.通常假定ε服从正态分布 N ( 0 , σ 2 ) N(0,σ2) N(0,σ2)

在这里插入图片描述
其中 σ \sigma σ为未知参数.称
y = a + b x + ε , ε ~ N ( 0 , σ 2 ) ( 1 ) y = a + b x +ε, ε ~N(0,σ2 )\qquad \qquad(1) y=a+bx+ε,εN(0,σ2)(1)
为一元线性回归模型. 由(1)得 E(y)=a+bx用E(y)作为y 的估计 y ^ \hat{y} y^
y ^ = a + b x ( 2 ) \hat{y}=a+bx \qquad\qquad (2) y^=a+bx(2)
称(2)为 y 关于 x 的一元线性回归方程 .

回归分析的任务是利用n组独立观察数据 ( x 1 , y 1 ) , … , ( x n , y n ) (x_1,y_1),…,(x_n,y_n) (x1,y1),,(xn,yn)来估计a和b, 以估计值 和 a ^ 和 b ^ \hat{a}和 \hat{b} a^b^分别代替(2)式中的a和b, 得回归方程

y ^ = a ^ + b ^ x \hat{y}=\hat{a}+\hat{b}x \qquad\qquad y^=a^+b^x

由于该方程的建立依赖于通过观察或试验取得的数据, 故又称其为经验回归方程或经验公式.
a ^ 和 b ^ \hat{a}和\hat{b} a^b^称为未知参数 a,b 的回归系数

最小二乘法估计

估计原则:
寻找一个使上述平方和达到最小的,作为这个物体重量的估计值, 这种方法称为最小二乘法
我们可以用真实值减去预测值的平方的和来检验回归的好坏,即
∑ i = 1 n ( y i − y ^ ) 2 \sum_{i=1}^{n}(y_i-\hat{y})^2 i=1n(yiy^)2
该值越小说明回归得到的效果越好。

依照最小二乘法的思想,提出目标量Q(a,b)
Q ( a , b ) = ∑ i = 1 n ( y i − ( a + b x ) ) 2 Q(a,b)=\sum_{i=1}^{n}(y_i-(a+bx))^2 Q(a,b)=i=1n(yi(a+bx))2
因此我们需要设法求出a 、b的估计值 a ^ 、 b ^ \hat{a}、\hat{b} a^b^,使偏差平方和Q(a,b)达到最小.
在这里插入图片描述
由此得到的回归直线 y ^ = a ^ + b ^ x \hat{y}=\hat{a}+\hat{b}x y^=a^+b^x 是在所有直线中偏差平方和Q(a,b)最小的一条直线.

由于 y ^ = a ^ + b ^ x \hat{y}=\hat{a}+\hat{b}x y^=a^+b^x 是从观察值得到的回归方程,它会随观察结果的不同而改变,并且它只反映了由 x 的变化引起的 y 的变化,并没有包含误差项.另外,能否给出未知参数 σ 2 \sigma^2 σ2的估计,由此引出以下问题:

  • σ 2 \sigma^2 σ2的点估计是什么?
  • 回归方程是否有意义? 即自变量 x 的变化是否真的对因变量 y 有影响? 因此有必要对回归效果作出检验.
  • 如果方程真有意义,用它预测y时,预测值与真值的偏差能否估计?

σ 2 \sigma^2 σ2的无偏估计

y ^ = a ^ + b ^ x i 称 为 y i − y ^ i 为 x i 处 的 残 差 称 Q e = ∑ i = 1 n ( y i − y ^ ) 2 = ∑ i = 1 n ( y i − a ^ − b ^ x i ) 2 \hat{y}=\hat{a}+\hat{b}x_i 称 为 y_i-\hat{y}_i为x_i 处的残差称Q_e=\sum_{i=1}^{n}(y_i-\hat{y})^2=\sum_{i=1}^{n}(y_i-\hat{a}-\hat{b}x_i)^2 y^=a^+b^xiyiy^ixiQe=i=1n(yiy^)2=i=1n(yia^b^xi)2为残差平方和
在这里插入图片描述
Q e Q_e Qe反映了除 x 外其它因素对 y 的影响, 这些因素没有反映在自变量x中, 它们可作为随机因素看待

线性假设的显著性检验(T 检验法)

在这里插入图片描述

r检验法

在这里插入图片描述

线性回归的方差分析(F 检验法)

在这里插入图片描述

利用回归方程进行预报(预测)

例1 为研究某一化学反应过程温度对产品得率的影响,测得观测数据如下:在这里插入图片描述

求Y关于x的回归方程u(x),并进行有关检验和预报

clc,clear;
x1=[100:10:190]';
y=[45 51 54 61 66 70 74 78 85 89]';
figure(1)
plot(x1,y,'*')%画数据散点图
hold on
x=[ones(10,1),x1];
[b,bint,r,rint,stats]=regress(y,x);
x2=[100:190];
plot(x2,b(1)+b(2)*x2)%画拟合图
figure(2)
rcoplot(r,rint)%画残差图
b,bint,stats

b =

   -2.7394
    0.4830


bint =

   -6.3056    0.8268
    0.4589    0.5072


stats =

   1.0e+03 *

    0.0010    2.1316    0.0000    0.0009

b - 多元线性回归的系数估计值
bint - 系数估计值的置信边界下限和置信边界上限
r - 残差
rint - 用于诊断离群值的区间
stats - 模型统计量,包括R^2 统计量、F 统计量及其 p 值,以及误差方差的估计值( σ 2 \sigma^2 σ2无偏估计)。
在这里插入图片描述

  • R平方决定系数,反应因bai变量的全部变异能通过回归关系被自变量du解释的比例。介于0~1之间,越接近1,回归拟合效果越好,一般认为超过0.8的模型拟合优度比较高。
  • F统计量是指在零假设成立的情况下,符合F分布的统计量。 F = ( R 2 / q ) / ( ( 1 − R 2 ) / ( n − k − 1 ) ) F=(R^2/q)/((1-R^2)/(n-k-1)) F=R2/q/1R2/(nk1))所以 R 2 = 0 时 , F = 0 ; R 2 = 1 时 , F = ∞ 。 R^2=0时,F=0;R^2=1时,F=\infin。 R2=0F=0;R2=1F=

可以看到,线性回归方程y=a+bx的a=-2.7394,b= 0.4830
a的置信区间为 [-6.3056, 0.8268]
b的置信区间为[ 0.4589 , 0.5072]
R 2 = 1 R^2=1 R2=1,p=0,说明回归效果显著

在这里插入图片描述
在这里插入图片描述

多元多项回归

应用条件

  1. 线性趋势:Y与Xi间具有线性关系
  2. 独立性:应变量Y的取值相互独立
  3. 正态性:对任意一组自变量取值,因变量Y服从正态分布
  4. 方差齐性:对任意一组自变量取值,因变量y的方差相同

后两个条件等价于:残差ε服从均数为0、方差为σ2的正态分布

一元多项式回归

在这里插入图片描述
在这里插入图片描述

在这里插入图片描述

例2

观测物体降落的距离s与时间t的关系,得到数据如下表,求s关于t的回归方程 s ^ = a + b t + c t 2 \hat{s}=a+bt+ct^2 s^=a+bt+ct2
在这里插入图片描述
方法一:直接作二次多项式回归(相当于拟合):

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];
plot(t,s,'*')
hold on
[p,S]=polyfit(t,s,2)
plot(t,p(1)*t.^2+p(2)*t+p(3))

在这里插入图片描述

方法二:化为多元线性回归

clc,clear
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);
plot(t,s,'*')
hold on
plot(t,b(3)*t.^2+b(2)*t+b(1))
b,stats
figure(2)
rcoplot(r,rint)

在这里插入图片描述

逐步回归分析

“最优”的回归方程就是包含所有对Y有影响的变量, 而不包含对Y影响不显著的变量回归方程.
选择“最优”的回归方程有以下几种方法:
(1)从所有可能的因子(变量)组合的回归方程中选择最优者;
(2)从包含全部变量的回归方程中逐次剔除不显著因子;
(3)从一个变量开始,把变量逐个引入方程;
(4)“有进有出”的逐步回归分析.

 以第四种方法,即逐步回归分析法在筛选变量方面较为理想.

逐步回归分析法的思想:

  • 从一个自变量开始,视自变量Y对作用的显著程度,从大到小地依次逐个引入回归方程.
  • 当引入的自变量由于后面变量的引入而变得不显著时,要将其剔除掉.
  • 引入一个自变量或从回归方程中剔除一个自变量,为逐步回归的一步.
  • 对于每一步都要进行Y值检验,以确保每次引入新的显著性变量前回归方程中只包含对Y作用显著的变量
  • 这个过程反复进行,直至既无不显著的变量从回归方程中剔除,又无显著变量可引入回归方程时为止

MATLAB实现

在这里插入图片描述
运行stepwise命令时产生三个图形窗口:Stepwise Plot,Stepwise Table,Stepwise History.
在Stepwise Plot窗口,显示出各项的回归系数及其置信区间.
Stepwise Table 窗口中列出了一个统计表,包括回归系数及其置信区间,以及模型的统计量剩余标准差(RMSE)、相关系数(R-square)、F值、与F对应的概率P.

例题

水泥凝固时放出的热量y与水泥中4种化学成分x1、x2、x3、 x4有关,今测得一组数据如下,试用逐步回归法确定一个 线性模型.
在这里插入图片描述

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];
stepwise(x,y)

运行结果如图
在这里插入图片描述
红色表示变量移出回归方程,绿色表示变量移入回归方程。当我们想看某变量与回归方程的关系时时可以点击使其变绿,其他变红即可。我们现在的目标是找rmse最小时,R方接近1时应该含有的变量,直接点击全部步骤就可以出来.

在这里插入图片描述
可以看到变量x1和x2的显著性最好,
对变量y和x1、x2作线性回归:

X=[ones(13,1) x1 x2];
b=regress(y,X)

b =
   52.5773
    1.4683
    0.6623

故最终模型为 y = 52.5773 + 1.4683 x 1 + 0.6623 x 2 y=52.5773+1.4683x_1+0.6623x_2 y=52.5773+1.4683x1+0.6623x2

多元线性回归的应用

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
提示可能存在多重共线性的情况:

  • 整个模型的检验结果为P<α,但各自变量的偏回归系数的检验结果P>α。
  • 专业上认为应该有统计学意义的自变量检验结果却无统计学意义。
  • 自变量的偏回归系数取值大小甚至符号明显与实际情况相违背,难以解释。
  • 增加或删除一个自变量或一条记录,自变量回归系数发生较大变化。

MATLAB实现

在这里插入图片描述

例3

设某商品的需求量与消费者的平均收入、商品价格的统计数据如下,建立回归模型,预测平均收入为1000、价格为6时 的商品需求量. 在这里插入图片描述
在这里插入图片描述

clc,clear
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')

方法一:直接用多元二项式回归

这里调用了多元线性回归的工具箱
在这里插入图片描述
代码中x1代表输入,x2代表输出,所以x1输入1000,x2输入6,得到需求量为 88.4791 ± 16.3387 88.4791\pm16.3387 88.4791±16.3387.
点击导出,可以将均方根误差也就是标准差rmse、残差residuals和beta导入到工作区查看
可以得到
beta =[110.5313, 0.1464, -26.5709, -0.0001, 1.8475]
rmse=4.5362
在这里插入图片描述
方法二: 化为多元线性回归

clc,clear
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=[ones(10,1) x1' x2' (x1.^2)' (x2.^2)'];
[b,bint,r,rint,stats]=regress(y,X);
b,stats


b =

  110.5313
    0.1464
  -26.5709
   -0.0001
    1.8475
stats =

    0.9702   40.6656    0.0005   20.5771

两个的结果相同,R方为0.97,说明置信度很高

非线性回归

在这里插入图片描述
通常选择的六类曲线如下:
在这里插入图片描述

观测物体降落的距离s与时间t的关系,得到数据如下表,求s关于t的回归方程 s ^ = a + b t + c t 2 \hat{s}=a+bt+ct^2 s^=a+bt+ct2
在这里插入图片描述
在这里插入图片描述
2.

function yhat=volum(beta,x)
     yhat=beta(1)*exp(beta(2)./x);
end
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]';
[beta,r,J]=nlinfit(x',y','volum',beta0);
[YY,delta]=nlpredci('volum',x',beta,r ,J);
plot(x,y,'k+',x,YY,'r')
beta

beta =
   11.6037
   -1.0641

在这里插入图片描述
结果为
在这里插入图片描述

练习题

在这里插入图片描述

clc,clear
x0=[20:5:65]';
y=[13.2 15.1 16.4 17.1 17.9 18.7 19.6 21.2 22.5 24.3]';
x=[ones(10,1),x0];
[b,bint,r,rint,stats]=regress(y,x);
hold on
plot(x0,y,'*',x0,b(1)+b(2)*x0)
figure(2)
rcoplot(r,rint)
b,bint,stats

在这里插入图片描述
在这里插入图片描述

clc,clear
x=[0:2:20];
y=[0.6 2 4.4 7.5 11.8 17.1 23.3 31.2 39.6 49.7 61.7];
%拟合的方法
a=polyfit(x,y,2);
y1=a(3)+a(2)*x+a(1)*x.^2;
%回归的方法
T=[ones(11,1) x' (x.^2)'];
[b,bint,r,rint,stats]=regress(y',T);
y2=b(1)+b(2)*x+b(3)*x.^2;
plot(x,y,'*',x,y1,'-.',x,y2,'ro')

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

function y=model(v,x)
b1=v(1);
b2=v(2);
b3=v(3);
b4=v(4);
b5=v(5);
y=(b1.*x(:,2)-(1/b5).*x(:,3)).*((1+b2.*x(:,1)+b3.*x(:,2)+b4.*x(:,3))).^(-1);
end
clc,clear
x1=[470 285 470 470 470 100 100 470 100 100 100 285 285];
x2=[300 80 300 80 80 190 80 190  300 300 80 300 190];
x3=[10 10 120 120 10 10 65 65 54 120 120 10 120];
x=[x1' x2' x3'];
y=[8.55 3.76 4.82 0.02 2.75 14.39 2.54 4.35 13.0 8.5 0.05 11.32 3.13]';
v=[1 0.05 0.02 0.1 2]';
[beta,r,j]=nlinfit(x,y,'model',v);
beta

beta =
    1.2046
    0.0605
    0.0383
    0.1084
    1.2399

在这里插入图片描述

function y=ln(beta,x)
y=beta(1)+beta(2).*log(x);
end
clc,clear
x=[2 3 4 5 7 9 12 14 17 21 28 56];
y=[35 42 47 53 59 65 68 73 76 82 86 99];
beta0=[12 2]';
[beta,r,J]=nlinfit(x',y','ln',beta0);
beta
plot(x,y,'o',x,beta(1)+beta(2).*log(x))

在这里插入图片描述

©️2020 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页