关闭

回归分析

标签: 统计学R语言
1651人阅读 评论(2) 收藏 举报
分类:

回归分析

1 一元线性回归


    1. 合金强度和碳含量

      由专业知识知道, 合金的强度 y 与合金中碳的含量有关, 为了解他们之间的关系, 从生产 中收集了一批数据, (xi,yi),i=1,,n 具体数据如下

      序号 碳含量x 强度 序号 碳含量x 强度y
      1 0.10 42.0 7 0.16 49.0
      2 0.11 43.5 8 0.17 53.0
      3 0.12 45.0 9 0.18 50.0
      4 0.13 45.5 10 0.20 55.0
      5 0.14 45.0 11 0.21 55.0
      6 0.15 47.5 12 0.23 60.0
               
  1. 数据读取与展示
    thl<-read.table("data/thl.csv",sep=",",header=TRUE)
    head(thl)
    
       thl   qd
    1 0.10 42.0
    2 0.11 43.5
    3 0.12 45.0
    4 0.13 45.5
    5 0.14 45.0
    6 0.15 47.5
    
  2. 散点图展示数据
    plot(thl,xlab="tanhanliang",ylab="Qiangdu")
    

    thl.png

  3. 一元回归分析模型
    Y=β0+β1X+ε
    • Y称为响应变量, X称为回归变量,
    • β1 为回归系数, β0 称为截距, 统称为回归 参数。
    • 参数 β1 表示如果自变量 x增加一个单位, 则响应变量 Y 相应增加 β1 个单位。
  4. 样本数据

    若 (xi,yi),i=1,,n 为 (X,Y) 的一组观测值, 则一元线性回归模型可以表示为

    yi=β0+β1xi+εi,i=1,,n

    其中误差相互独立, 并满足 E(εi)=0,var(εi)=σ2,i=1,,n

  5. 回归参数的估计-最小二乘法

    目标函数为(误差平方和)

    Q(β0,β1)=i=1n(yiβ0β1xi)2

    使上述目标函数达到最小的 β0,β1 称为参数的最小二乘估计 β^0,β^1 ,即

    Q(β^0,β^1)=minβ0β1Q(β0,β1)
  6. 求解最小二乘问题

    关于 β0,β1 分别求导, 并令之为0, 得方程组如下

    β0Q(β0,β1)β1Q(β0,β1)=2i=1n(yiβ0β1xi)=0=2i=1nxi(yiβ0β1xi)=0

    上述方程组可以写作如下形式, 称为正规方程组

    2i=1n(1xi)(yiβ0β1xi)=(00)
  7. 最小二乘法所得参数的估计
    β^1=i=1n(xix¯)(yiy¯)i=1n(xix¯)2=SxySxx,β^0=y¯β^1x¯
    • 其中
    x¯=1ni=1nxi,y¯=1ni=1nyi,Sxx=i=1n(xix¯)2,Syy=i=1n(yiy¯)2,Sxy=i=1n(xix¯)(yiy¯)
  8. 上述诸变量的计算
    • 在进行最小二乘法进行手动计算时,比较容易得出的是
      n,i=1nxi,i=1nyi,i=1nx2i,i=1nxiyi
    • 下面给出上述估计量的计算
      x¯SxxSxy=1ni=1nxi,y¯=1ni=1ni=1nyi=i=1n(xix¯)(xix¯)=i=1n(xix¯)xix¯i=1n(xix¯)=i=1nx2inx¯2=i=1n(xix¯)(yiy¯)=i=1n(xix¯)yi=i=1nxiyinx¯y¯
  9. 误差方差的估计

    注意到 σ2 是未知的, 通过拟合之后的残差可以得出 σ2 的估计为

    σ^2=1n2i=1ne^2i

    作为参数 σ2 的估计,

    • 其中 e^i=yiβ^0β^1xi 为残差
    • 可以证明 σ^2 为 σ2 的无偏估计。
  10. 可以得出参数估计的均值和方差
    E(β^0)=β0,E(β^1)=β1

    说明估计量是无偏的,通过计算可得截距和斜率的估计的方差为

    Var(β^0)=σ2(1n+x¯2Sxx),Var(β^1)=σ2Sxx
  11. 参数估计量的分布

    若给定误差服从正态分布, 可以得到

    β^1N(β1,σ2Sxx),(β^1β1)Sxxσ2N(0,1)

    另外有:

    (n2)σ^2σ2χ2(n2)

    可以得出

    (β^1β1)Sxxσ^2t(n2)

    根据这个分布可以进行回归方程的显著性检验

  12. 回归方程的显著性检验
    • 从统计上讲, β1 是 E(Y) 随 X 线性变化的变化率, 若 β1=0 , 回归方程没有意义
    • 回归方程的显著性检验需要检验如下假设:
      H0:β1=0,H1:β10
    • 原假设成立时( β1=0 )可以采用t检验统计量对模型显著性进行检验
      t=β^1Sxxσ^t(n2)

      或 F 检验统计量

    F=t2=β1^2Sxxσ^2F(1,n2)
  13. 教材例子
    1. 炼钢厂出钢时采用的盛钢水的钢包

      在使用过程中, 由于钢液及炉渣对包衬耐火材料的侵蚀, 使其容积不断增大, 希望找出使用 次数与增大的容积之间的关系

  14. 散点图
      x=c(2:16);
      y=c(6.42,8.20,9.58,9.50,9.70,10,9.93,9.99,10.49,10.59,10.6,10.8,10.6,10.9,10.76)
    plot(x,y,xlab="使用次数",ylab="增大容积")
    

    gangye.png

  15. 建模

    根据数据特点, 首先进行非线性变换, 然后拟合如下模型

    lnY=lna+bX

    令 Y=ln(Y),X=X1,β0=ln(a) , 则可以化为线性模型

    Y=β0+β1X
    ynew=log(y);xnew=1/x; plot(xnew,ynew)
    

    gangye2.png

  16. 线性回归过程
    (l<-lm(ynew~xnew))
    
    Call:
    lm(formula = ynew ~ xnew)
    
    Coefficients:
    (Intercept)         xnew
           2.46        -1.11
    
  17. 线性回归详细结果
    summary(l)
    
    Call:
    lm(formula = ynew ~ xnew)
    
    Residuals:
         Min       1Q   Median       3Q      Max
    -0.04303 -0.01506  0.00310  0.00611  0.07956
    
    Coefficients:
                Estimate Std. Error t value Pr(>|t|)
    (Intercept)   2.4578     0.0126   195.2  < 2e-16 ***
    xnew         -1.1107     0.0638   -17.4  2.2e-10 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    Residual standard error: 0.029 on 13 degrees of freedom
    Multiple R-squared: 0.959,	Adjusted R-squared: 0.956
    F-statistic:  303 on 1 and 13 DF,  p-value: 2.17e-10
    
  18. 拟合后的图形
    plot(xnew,ynew)
    abline(l,col=2)
    

    gangye4.png

  19. matlab 实现上述过程
    x=[2:16]';
    y=[6.42;8.20;9.58;9.50;9.70;10;9.93;9.99;
       10.49;10.59;10.6;10.8;10.6;10.9;10.76];
    ynew=log(y);xnew=[ones(15,1),1./x];
    [b,bint,~,~,stats]=regress(ynew,xnew)
    
    b =
    2.4578
    -1.1107
    bint =
    2.4306    2.4850
    -1.2485   -0.9729
    stats =
    0.9589  303.1896    0.0000    0.0008
    

2 (附录)一元线性回归的求解过程

  1. 回归参数的最小二乘法估计

    目标函数为(误差平方和)

    Q(β0,β1)=i=1n(yiβ0β1xi)2

    关于 β0,β1 分别求导, 并令之为0, 得方程组如下

    β0Q(β0,β1)β1Q(β0,β1)=2i=1n(yiβ0β1xi)=0=2i=1nxi(yiβ0β1xi)=0

    上述方程组可以写作如下形式, 称为正规方程组

    2i=1n(1xi)(yiβ0β1xi)=(00)
  2. 正规方程组解的存在性

    将上述方程组进行化简

    nβ0+(i=1nxi)β1=i=1nyi,(i=1nxi)β0+(i=1nx2i)β1=i=1nxiyi

    称为正规方程组, 由于 xi 不完全相同, 其系数行列式

    ni=1nxii=1nxii=1nx2i=ni=1n(xix)20

    上述方程有唯一解

  3. 正规方程组
    2i=1n(1xi)(yiβ0β1xi)=(00)

    将上述方程展开,写作

    2111x1x2xnTy1β0β1x1y2β0β1x2ynβ0β1xn=(00)
  4. 方程组的矩阵表示:
    • Y=y1y2yn,X=111x1x2xn,β=(β0β1)
    • 则上述方程可以简化为
      2XT(YXβ)=0
  5. 方程组的另一种解释(向量求导数)
    Q(β)=Q(β0,β1)=(YXβ)T(YXβ)

    则关于 β 求偏导数, 令之为零, 有

    βQ(β)=2XT(YXβ)=0

    求解

    XTXβ=XTY
  6. 回归参数的最小二乘估计

    解之可得:

    β^=(XTX)1XTY=nnx¯nx¯i=1nx2i1i=1nyii=1nxiyi
  7. 回归参数估计的性质

    注意到

    β^=(XTX)1XTY,Y=Xβ+ε

    β^=β+(XTX)1XTε

    由此可得估计量的均值和方差为

    E(β^)=β,Var(β^)=σ2(XTX)1
  8. 截距和斜率估计的方差

    注意到:

    (XTX)1=nnx¯nx¯i=1nx2i1=1n+x¯2Sxxx¯Sxxx¯Sxx1(xix¯)2

    可得 β^0,β^1 的方差分别为

    Var(β^0)=σ2(1n+x¯2Sxx),Var(β^1)=σ2Sxx

3 多元回归与逐步回归

  1. 多元回归分析的数学模型

    设自变量为 x1,x2,,xp 目标变量为 Y , 满足

    y=β0+β1x1++βpxp+ε

    设从上述模型得到 n 个相互独立的观测,即有

    y1=β0+β1x11+β2x12++βpx1p+ε1y2=β0+β1x21+β2x22++βpx2p+ε2