基于R的高级统计（全局优化）

作者: Shaoqian_Ma | 来源:发表于2020-04-23 11:42 被阅读0次

General Optimization

全局优化，给定f和初始的x，希望找到f的最小值

在k维空间中选择一个方向pn；
选择在pn方向的一个步长α，使得：
$min_{α>0}f(x_n+αp_n)$
从而更新我们对x的估计：
$x_{n+1}=x_n+α_np_n$

解决这样的问题需要让每次的步骤的效率尽可能高，因为k可能是个很大的数，有很多未知的参数需要优化（高维参数空间），所以值得考虑的，一方面减少内存消耗，一方面减少计算时间

Steepest Descent

直译过来是最陡下降算法，或者叫梯度下降法，是一种一阶最优化算法。用这个方法来解决从什么方向下降，以及每次下降的步长的问题。

直观地讲，我们当然希望下降的方向是“最陡“，f变化最快的那个，反映在等高线上就是与xn所在切线垂直的方向。f的变化快慢用导数来衡量，因此有：
$x_{n+1}=x_n−α_nf′(x_n)$
但是每次都选择最陡的方向在某些情况下可能会有问题：比如说参数之间有相关性，比如下图的椭圆等高线

steepestdescent2.png

这样在到达最小值点的过程中就需要很多步，反而是在走弯路。另外梯度下降需要的步数与初始值选取关系很大。

Example: Multivariate Normal

以多元正态分布为例，计算负对数似然：

set.seed(2017-08-10)
mu <- c(1, 2)#二元正态分布的均值，也就是需要在等高线图中找到的那个点
S <- rbind(c(1, .9), c(.9, 1))#协方差矩阵
x <- MASS::mvrnorm(500, mu, S)#创建数据集，#生成多元正态数据，使用MASS 包中的 mvrnorm() 函数，其格式为mvrnorm(n, mean, sigma)，
#其中 n 是你想要的样本大小， mean 为均值向量，而 sigma 是方差—协方差矩阵（或相关矩阵）
nloglike <- function(mu1, mu2) {
        dmv <- mvtnorm::dmvnorm(x, c(mu1, mu2), S, log = TRUE)#对于x每一个点（x1,x2)，计算概率密度的负对数似然（因为是点所以用密度近似），对数求和sum就是概率相乘。
        -sum(dmv)
}
nloglike <- Vectorize(nloglike, c("mu1", "mu2"))#方便可视化
nx <- 40
ny <- 40
xg <- seq(-5, 5, len = nx)
yg <- seq(-5, 6, len = ny)
g <- expand.grid(xg, yg)#生成网格点，相当于固定一个x1轴，另一个轴x2生成很多点
nLL <- nloglike(g[, 1], g[, 2])#对数据集x，按均值向量里每一个(mu1,mu2)返回对应的负对数似然值作为空间坐标里的z轴
z <- matrix(nLL, nx, ny)#(40*40)格点的矩阵，nll相当于立体空间里的高度值
par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), #expression方便在坐标上打印表达式，比如mu1
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)

example_1.png

最值点在（1,2），假设我们从（-5，-2）点起始开始梯度下降并记录下降路径：

library(dplyr, warn.conflicts = FALSE)
norm <- function(x) x / sqrt(sum(x^2))
Sinv <- solve(S)  ## I know I said not to do this!实际运算中尽量不要去求矩阵的逆
step1 <- function(mu, alpha = 1) {
        D <- sweep(x, 2, mu, "-")#对数组或者矩阵进行运算。 MARGIN=1表示行，2表示列，函数默认是减法
#-Usage
#sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...)
#对数组的某一个或某几个维度减去（或FUN指定的操作）STATS
        score <- colSums(D) %>% norm
        mu + alpha * drop(Sinv %*% score)
}
steep <- function(mu, n = 10, ...) {
        results <- vector("list", length = n)
        for(i in seq_len(n)) {
                results[[i]] <- step1(mu, ...)
                mu <- results[[i]]
        }
        results
}
m <- do.call("rbind", steep(c(-5, -2), 8))
m <- rbind(c(-5, -2), m)

par(mar = c(4.5, 4.5, 1, 1))
contour(xg, yg, z, nlevels = 40, xlab = expression(mu[1]), 
        ylab = expression(mu[2]))
abline(h = 0, v = 0, lty = 2)
points(m, pch = 20, type = "b")

example_2.png

可以看到实际路径是曲折下降的

The Newton Direction

对f用多项式近似：
$f(x_n+p)≈f(x_n)+p′f′(x_n)+\frac 12p′f′′(x_n)p$
把右侧看成是对p的函数，求其极小值点，得到：
$p_n=f′′(x_n)^{−1}[−f′(x_n)]$
因此有迭代公式：
$x_{n+1}=x_n−f′′(x_n)^{−1}f′(x_n)$
这里右侧的乘积是：Hessian矩阵的逆乘以f对每个元素求偏导。

每一次迭代，都相当于把对f进行优化的问题转化成用泰勒展开的g来近似f，求g的最值。

直观地来看：

curve(-dnorm(x), -2, 3, lwd = 2, ylim = c(-0.55, .1))
xn <- -1.2
abline(v = xn, lty = 2)
axis(3, xn, expression(x[n]))
g <- function(x) {#x-xn就是p，对正态概率密度函数求导有φ'(x)=(-x)*φ(x)
        -dnorm(xn) + (x-xn) * xn * dnorm(xn) - 0.5 * (x-xn)^2 * (dnorm(xn) - xn * (xn * dnorm(xn)))
}
curve(g, -2, 3, add = TRUE, col = 4)
op <- optimize(g, c(0, 3))#优化函数，(0,3)是x取值区间，
abline(v = op$minimum, lty = 2)
axis(3, op$minimum, expression(x[n+1]))

Newtom_direction.png

这个图黑色的线是原密度函数，蓝色的线是用牛顿近似的。可以看出，仅一次步长就使得更新后的x已经远离了原来的最小值点。牛顿法的典型特点就是它并不能保证每一步都使得你距离希望的那个点更近，结果不稳定。但这个方法确实速度是最快的（因为太快所以很容易越过真正的极值点）。后面我们会讨论的EM算法与此相反，速度很慢但保证每一步都是在提升。

但下图我们发现再做第二次近似时，又往极值点更近了一些：

curve(-dnorm(x), -2, 3, lwd = 2, ylim = c(-0.55, .1))
xn <- -1.2
op <- optimize(g, c(0, 3))
abline(v = op$minimum, lty = 2)
axis(3, op$minimum, expression(x[n+1]))

xn <- op$minimum
curve(g, -2, 3, add = TRUE, col = 4)
op <- optimize(g, c(0, 3))
abline(v = op$minimum, lty = 2)
axis(3, op$minimum, expression(x[n+2]))

Newton_direction2.png

Generalized Linear Models

广义线性模型：是标准线性模型在非正态响应变量上应用的扩展。响应变量服从指数族分布，通常写作：
$y_i∼p(y_i∣μ_i)$
其中p为指数族分布，yi的期望为μi

应用GLM的设计假设：https://blog.csdn.net/TRillionZxY1/article/details/77140539 理解自然参数η和μ的关系

假设函数和模型比较像，一般就可以理解成还没加工好的模型，参考：https://blog.csdn.net/weixin_40166430/article/details/81211822 有了假设函数就可以构造损失函数
$g(μ_i)=x′_iβ$
上式中g是一个非线性函数，是一种link funciton，也就是η=η(μ)的反函数
$Var(y_i)=V(μ)$
V是一个已知的方差函数（variance function)

与标准线性模型不同，非线性对参数β的最大似然估计不是closed form，也就是不能显式地用最小二乘法代入变量求得参数的值，而是必须通过迭代的方法来获取估计。

经典的方法是Fisher score algorithm，去线性近似g，写成下式：其实也是泰勒展开
$g(y_i)≈g(μ_i)+(y_i−μ_i)g′(μ_i)$
其原理如下：了解一下即可

GLM_Fisher.jpg

举例：Poisson回归

响应变量y满足：
$y_i∼Poisson(μ_i)$

$g(μ)=logμ_i=x′_iβ$

这里的g就是link function，简单理解为：Y和X本来成指数关系则给Y取对数得到log Y和X成线性关系，才可以建立对应的线形模型，是为glm，其中的log就是link function。建立η=xβ和期望μ之间的关系

Newton’s Method in R

在R里面可以用 nlm() 函数在给定一个初始值向量的情况下实现牛顿法优化目标函数。logistic模型即是在伯努利分布和广义线性模型的假设下推导而来，其参数φ(p)就是sigmoid函数，属于广义线性模型。这里以logistic回归为例：
$y_i∼Bernoulli(p_i)$
对数线性模型为：
$log\frac {p_i}{1−p_i}=β_0+x_iβ_1$
目标是利用最大似然法估计β的值。
$logL(β)=∑_{i=1}^ny_i(β_0+x_iβ_1)−log(1+e^{(β_0+x_iβ_1)})$
为了利用牛顿法使对数似然取最小，对β求导：
$ℓ′(β)=∑_{i=1}^nℓ′_i(β)$

$ℓ′′(β)=∑_{i=1}^nℓ′′_i(β)$

这样接下来就可以应用牛顿法求参数最优解了。在R里可以使用deriv()函数计算导数，

Example: Trends in p-values Over Time

使用到的包是：

remotes::install_github("jtleek/tidypvals")

目的是看p-value与年份的回归关系，把响应变量p-value以0.05为界限划分成两类，预测变量x是year，减去2000是因为我们希望以2000年为起点

library(tidypvals)
library(dplyr)
jager <- mutate(tidypvals::jager2014, 
                pvalue = as.numeric(as.character(pvalue)),
                y = ifelse(pvalue < 0.05 
                           | (pvalue == 0.05 & operator == "lessthan"), 
                           1, 0),
                x = year - 2000) %>%
        tbl_df

接下来计算负对数似然对β0和β1的梯度和Hessian矩阵。在deriv函数中，声明β0和β1是变量，function.arg = TRUE是为了返回一个以β为参数的函数

nll_one <- deriv(~ -(y * (b0 + x * b1) - log(1 + exp(b0 + b1 * x))),
             c("b0", "b1"), function.arg = TRUE, hessian = TRUE)

nll_one
function (b0, b1) 
{
    .expr6 <- exp(b0 + b1 * x)
    .expr7 <- 1 + .expr6
    .expr11 <- .expr6/.expr7
    .expr15 <- .expr7^2
    .expr18 <- .expr6 * x
    .expr19 <- .expr18/.expr7
    .value <- -(y * (b0 + x * b1) - log(.expr7))
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("b0", 
        "b1")))
    .hessian <- array(0, c(length(.value), 2L, 2L), list(NULL, 
        c("b0", "b1"), c("b0", "b1")))
    .grad[, "b0"] <- -(y - .expr11)
    .hessian[, "b0", "b0"] <- .expr11 - .expr6 * .expr6/.expr15
    .hessian[, "b0", "b1"] <- .hessian[, "b1", "b0"] <- .expr19 - 
        .expr6 * .expr18/.expr15
    .grad[, "b1"] <- -(y * x - .expr19)
    .hessian[, "b1", "b1"] <- .expr18 * x/.expr7 - .expr18 * 
        .expr18/.expr15
    attr(.value, "gradient") <- .grad
    attr(.value, "hessian") <- .hessian
    .value
}

这样，nll_one函数就可以通过β传参，计算其对应的负对数似然。而且上面函数的属性里包括了梯度和Hessian矩阵。以我们自己的数据集为例：

x <- jager$x
y <- jager$y
str(nll_one(0, 0))

num [1:15653] 0.693 0.693 0.693 0.693 0.693 ...
 - attr(*, "gradient")= num [1:15653, 1:2] -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : NULL
  .. ..$ : chr [1:2] "b0" "b1"
 - attr(*, "hessian")= num [1:15653, 1:2, 1:2] 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
  ..- attr(*, "dimnames")=List of 3
  .. ..$ : NULL
  .. ..$ : chr [1:2] "b0" "b1"
  .. ..$ : chr [1:2] "b0" "b1"

上面打印的是nll_one函数对于x,y的每一个点计算了负对数似然，但是没有加和。所以需要额外求加和的话可以自己写个函数：

nll <- function(b) {
        v <- nll_one(b[1], b[2])
        f <- sum(v)                                     ## log-likelihood
        gr <- colSums(attr(v, "gradient"))              ## gradient vector
        hess <- apply(attr(v, "hessian"), c(2, 3), sum) ## Hessian matrix
        attributes(f) <- list(gradient = gr, 
                              hessian = hess)
        f
}

测试一下，假设β全零：

nll(c(0, 0))
[1] 10849.83
attr(,"gradient")
      b0       b1 
 -4586.5 -21854.5 
attr(,"hessian")
         b0        b1
b0  3913.25  19618.25
b1 19618.25 137733.75

以β0=0，β1=0为起始点，就可以用nlm函数对负对数似然求最小值：

res <- nlm(nll, c(0, 0))
res
$minimum
[1] 7956.976

$estimate
[1]  1.57032807 -0.04416515

$gradient
[1] -0.000001451746 -0.000002782241

$code
[1] 4

$iterations
[1] 100

code值为4表示迭代次数超过了上限，在这里也就是iterations=100，说明函数可能还没有收敛，不过我们可以设置迭代上限：

res <- nlm(nll, c(0, 0), iterlim = 1000)
res
$minimum
[1] 7956.976

$estimate
[1]  1.57032807 -0.04416515

$gradient
[1] -0.000001451746 -0.000002782241

$code
[1] 2

$iterations
[1] 260

code值为2表面迭代结果不错，260次迭代低于我们设置的上限1000.

实际上大多数的优化算法都会有个可选的设置来设定传入参数相对的大小，一般最后输出的参数的数量级是差不多的。当我们的函数中包含的参数大小差别较大，不在一个数量级时，对计算机而言可能会存在一些问题。在nlm函数中可以设置typsize参数，与参数长度等长，设置参数彼此间相对的大小：

res <- nlm(nll, c(0, 0), iterlim = 1000,
           typsize = c(1, 0.1))#表示β0大致是β1的10倍
res
$minimum
[1] 7956.976

$estimate
[1]  1.57032807 -0.04416515 

$gradient
[1] -0.000001451745 -0.000002782238

$code
[1] 1

$iterations
[1] 4

虽然结果是一样的，但是迭代次数明显减少了，iterations = 4。在实际应用中实现给定合适的参数大小往往可以使得优化算法表现得更好。

参考教材：Advanced Statistical Computing

基于R的高级统计（全局优化）

General Optimization

Steepest Descent

Example: Multivariate Normal

The Newton Direction

Generalized Linear Models

Newton’s Method in R

Example: Trends in p-values Over Time

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读