[关键字]: R, differentiation, calculas, algebraic computation
[出处]: http://rpubs.com/wush978/differentiation
[重点摘要]:
ps. 好读版请参考上面的link。
最近笔者有计算一个函数的Gradient和Hessian的需求，所以研究了一下R 中
进行这类计算的功能和套件。
以下我会介绍R 中进行数值微分和代数微分的功能。
# numDeriv
我第一个介绍的是数值微分套件：`numDeriv`。他的用法很简单：
## `grad`
`grad`函数可以拿来计算函数在指定点的gradient：
```{r}
library(numDeriv)
f <- function(x) x^2
grad(f, 2)
```
`grad`也可以计算向量函数的gradient：
```{r}
f <- function(x) x[1] * x[2]
grad(f, c(1,2))
```
## hessian
`hessian`可以计算函数在指定点的hessian
```{r}
hessian(f, 1:2)
```
由于是数值解的关系，所以可以看到二阶微分会算出很小的值。
# `deriv`
接下来我要介绍一个做代数微分的功能。
首先我们把刚刚的`f <- function(x) x^2`写成expression：
```{r}
f <- function(x) x^2
f.exp <- expression(x^2)
```
接着使用内建的`deriv`：
```{r}
deriv(f.exp, namevec="x")
```
R 就会回传出计算gradient的方式。将它复制贴上建立函数，并且修改内容
让它是vectorized：
```{r}
g <- function(x) {
.value <- x^2
.grad <- 2 * x
.grad
}
```
和`numDeriv`的结果做个比较：
```{r}
all.equal(g(2), grad(f, 2))
```
`deriv`的另一个参数`hessian`设为`TRUE`的时候，回传值会附加hessian的
算式。
# 范例
根据`r citet(bib["bowling2009logistic"])`，以下函数是standard
normal cdf的logistic approximation：
$$F(x) = \frac{1}{1+e^{-0.07056 x^3 - 1.5976 x}}$$
```{r}
F.0 <- function(x) 1/(1 + exp(-0.07056 * x^3 - 1.5976 * x))
curve(pnorm, -3, 3, lty = 1)
curve(F.0, add = TRUE, col = 2, lty = 2)
```
肉眼根本分不出来差异！
如果我们用它来取代某些Likelihood中的standard normal cdf(当资料有
censoring的时候)，在计算MLE的时候就会需要计算$F'(x)$和$F''(x)$。他
们分别是：
- $$F'(x) =
-{{\left(-0.21168\,x^2-1.5976\right)\,e^{-0.07056\,x^3-1.5976\,x}
}\over{\left(e^{-0.07056\,x^3-1.5976\,x}+1\right)^2}}$$
- $$F''(x) =
-{{\left(-0.21168\,x^2-1.5976\right)^2\,e^{-0.07056\,x^3-1.5976\,x
}
}\over{\left(e^{-0.07056\,x^3-1.5976\,x}+1\right)^2}} +
\frac{0.42336\,x\,
e^{-0.07056\,x^3-1.5976\,x}}{\left(e^{-0.07056\,x^3-1.5976\,x} +
1\right)^2} +
{{2\,\left(-0.21168\,x^2-1.5976\right)^2\,e^{-0.14112\,
x^3-3.1952\,x}}\over{\left(e^{-0.07056\,x^3-1.5976\,x}+1\right)^3}
}$$
好丑好难写。
这时候就可以利用`deriv`来建立gradient and hessian。先建立gradient：
```{r}
F.exp <- expression(1/(1 + exp(-0.07056 * x^3 - 1.5976 * x)))
deriv(F.exp, namevec="x")
```
稍作修改就可以得到：
```{r}
F.1 <- function(x) {
.expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
.expr7 <- 1 + .expr6
.value <- 1/.expr7
.expr6 * (0.07056 * (3 * x^2) + 1.5976)/.expr7^2
}
```
再参考`deriv`的hessian算式：
```{r}
deriv(F.exp, namevec="x", hessian=TRUE)
```
可以写出：
```{r}
F.2 <- function(x) {
.expr6 <- exp(-0.07056 * x^3 - 1.5976 * x)
.expr7 <- 1 + .expr6
.expr12 <- 0.07056 * (3 * x^2) + 1.5976
.expr13 <- .expr6 * .expr12
.expr14 <- .expr7^2
1/.expr7
(.expr6 * (0.07056 * (3 * (2 * x))) -
.expr13 * .expr12)/.expr14 + .expr13 * (2 * (.expr13 *
.expr7))/.expr14^2
}
```
和`numDeriv`比较一下：
```{r}
x <- rnorm(1)
grad(F.0, x)
F.1(x)
hessian(F.0, x)
F.2(x)
```
完全一致！
但是用`deriv`建立的函数效能绝对比较好：
```{r}
library(microbenchmark)
microbenchmark(grad(F.0, x), F.1(x))
microbenchmark(hessian(F.0, x), F.2(x))
```
以上内容给大家参考囉！
# Reference
```{r, results='asis', echo=FALSE}
bibliography()
```