美文网首页
PH525x series - Introduction to

PH525x series - Introduction to

作者: 3between7 | 来源:发表于2019-11-22 17:38 被阅读0次

    本章会对线性模型做一个大致的介绍,还是举例说明吧:

    例1:自由落体问题

    想象自己是16世纪的伽利略,正在研究自由落体问题,一名助理爬上了比萨斜塔扔下了一个球,同时还另有几名助理记录不同时间球所在的位置。现借助如今我们已知的自由落体公式(加上未知的测量错误)来模拟圆球的运动轨迹:

    set.seed(1)
    g <- 9.8 ##meters per second
    n <- 25
    tt <- seq(0,3.4,len=n) ##time in secs, note: we use tt because t is a base function
    d <- 56.67  - 0.5*g*tt^2 + rnorm(n,sd=1) ##meters
    mypar()
    plot(tt,d,ylab="Distance in meters",xlab="Time in seconds")
    
    Galileo.png

    虽然伽利略并不知道自由落体公式,但通过描点绘图依旧可以得到上图,从图中可以看出,这条轨迹形似抛物线,所以伽利略建立如下等式:

    Y_i = β_0 + β_1x_i + β_2x_i^2 + ε,i=1,...,n

    其中,Y_i代表圆球距离地面的高度,x_i代表下落时间,ε代表测量误差。因为这个方程式是已知统计量与未知参数的线性组合,所以被称为线性模型。

    例2:父与子身高问题

    现在想象自己是19世纪的高尔顿,收集到了很多对父子的身高数据,其身高分布形如:

    scatterPlot.png

    看图可知,儿子的身高与父亲的身高之间大致呈现出了一个线性正相关的趋势。本例中,可用如下模型描述父与子的身高数据:
    Y_i = β_0 + β_1x_i + ε,i=1,...,N

    这同样是x_iY_i的线性模型,其中,模型固定了父亲的身高,所以是小写。另外,单单是测量误差并不能完全解释ε中的变量,说明还有其他变量未被纳入到这个模型中,比如说,母亲的身高、遗传随机性以及环境因子等。

    例3:来自不同群体的随机样本问题

    我们使用小鼠体重数据(饲喂了两种不同饲料),每组随机选择12只小鼠。小鼠体重数据分布如下:

    mice.png

    为了研究不同群体间平均体重的差异,除了使用t检验外,还可以用线性模型:

    Y_i = β_0 + β_1x_i + ε_i

    其中,β_0是chow组的平均体重,β_1是两组体重均值之差,当第i只小鼠的饲料分别是是hf、chow时,x_i分别等于1和0,而ε_i是同一群体内小鼠间的差异。

    一般线性模型

    一般线性模型的方程式如下:

    Y_i = β_0 + β_1x_{i,1} + β_2x_{i,2} + ... + β_px_{i,p} + ε,i=1,...,n
    Y_i = β_0 + \sum_{j=1}^pβ_jx_{i,j} + ε,i=1,...,n

    矩阵代数提供了一种简洁的语言和数学框架,可以用任何符合上述框架的线性模型进行计算和推导。

    参数估计

    若想让上述线性模型有意义,需要去估计未知参数β的值。在一个例子中,我们想要描述一个物理过程,所以不能有未知参数;第二个例子里,我们想要探究父亲的身高在平均水平上可多大程度影响儿子的身高;而在最后一个例子中,我们想要高明白两个群体间的体重实际上是否有差异,也就是是否β_1≠0

    科学的一般做法是找到可以最小化拟合模型与实际数据之间的距离的值,下面这个表达式叫做最小二乘方程(LS)
    \sum_{i=1}^n\{Y_i - (β_0 + \sum_{j=1}^pβ_jx_{i,j})\}^2

    其实就是求ε的平方和啥时候最小,一旦我们发现了最小值,就可以将这个值叫做最小二乘估计(LSE),记为\hatβ,在估计过程中求最小二乘方程时得到的统计量被称为残差平方和(RSS)。由于这些统计量都依赖Y值,所以它们都属于随机变量。

    lm函数求LSE

    lm函数的作用就是进行线性拟合,仍拿第一个例子说明,经过绘图可以发现图形的轨迹类似抛物线,所以可以使用如下表达式进行拟合:
    Y_i = β_0 + β_1x_i + β_2x_i^2 + ε,i=1,...,n
    使用R语言的lm()函数便可实现这一拟合过程,返回结果中的coef就是LSE值:

    tt2 <-tt^2
    fit <- lm(y~tt+tt2)
    summary(fit)$coef
    
    ##               Estimate Std. Error    t value     Pr(>|t|)
    ## (Intercept) 57.1047803  0.4996845 114.281666 5.119823e-32
    ## tt          -0.4460393  0.6806757  -0.655289 5.190757e-01
    ## tt2         -4.7471698  0.1933701 -24.549662 1.767229e-17
    

    那么,lm函数是如何计算LSE的呢?为了回答这一问题,首先,让我们编写一个给定β向量计算RSS的函数:

    rss <- function(Beta0,Beta1,Beta2){
      r <- y - (Beta0+Beta1*tt+Beta2*tt^2)
      return(sum(r^2))
    }
    

    所以,给定任意三维向量,我们便可以计算RSS的值。接下来我们固定β_0β_1的值,将该其变为β_2的函数:

    Beta2s<- seq(-10,0,len=100)
    plot(Beta2s,sapply(Beta2s,rss,Beta0=55,Beta1=0),
         ylab="RSS",xlab="Beta2",type="l")
    ##Let's add another curve fixing another pair:
    Beta2s<- seq(-10,0,len=100)
    lines(Beta2s,sapply(Beta2s,rss,Beta0=65,Beta1=0),col=2)
    
    lm.png

    在这里使用试错法是行不通的,相反,我们可以使用微积分:取偏导数,设为0,然后求解。

    偏导数:在数学中,一个多变量的函数的偏导数,就是它关于其中一个变量的导数而保持其他变量恒定

    “取偏导数,设为0,然后求解。”这句话具体到自由落体这个例子中,我的理解就是:固定了其他变量之后,求曲线上导数为0,也就是斜率为0的那个点,该点对应的纵坐标便是一个个RSS,再从这一堆RSS中找到最小值。

    参考文章

    相关文章

      网友评论

          本文标题:PH525x series - Introduction to

          本文链接:https://www.haomeiwen.com/subject/kmtqwctx.html