作者Wush978 (拒看低質媒體)
看板R_Language
標題[分享] R 中計算函數的微分
時間Thu Jan 16 00:29:16 2014
[關鍵字]: 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()
```
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 1.34.138.85
1F：推 maninblue:推~ 01/16 11:11