deriv を利用して、R で記号微分してみます。
正規分布の対数尤度関数を記号微分してみよう
正規分布の対数尤度関数は、
\[
l(\mu, \s; x)=-\frac{\log(2\pi\s^2)}{2}-\frac{(x - \mu)^2}{2\s^2}
\] ですので、これを偏微分すると、
\[
\begin{align*}
\frac{\partial l(\mu, \s; x)}{\partial \mu} &= \frac{x-\mu}{\s^2} \\
\frac{\partial l(\mu, \s; x)}{\partial \s} &= -\frac{1}{\s}+\frac{(x-\mu)^2}{\s^3}
\end{align*}
\] となる予定です。
deriv を用いた数値微分では、まず被微分関数を expression、call、または formula の形で定義します。
fx <- expression( -log(2*pi*sd^2)/2 - (x - mean)^2/(2*sd^2) )
あとは、
を指定して実行するだけです。
dfx <- deriv( fx, c("mean", "sd"), function.arg = c("x", "mean", "sd") ) dfx
結果は以下のようになりました。
> dfx function (x, mean, sd) { .expr1 <- 2 * pi .expr2 <- sd^2 .expr3 <- .expr1 * .expr2 .expr7 <- x - mean .expr8 <- .expr7^2 .expr9 <- 2 * .expr2 .expr14 <- 2 * sd .value <- -log(.expr3)/2 - .expr8/.expr9 .grad <- array(0, c(length(.value), 2L), list(NULL, c("mean", "sd"))) .grad[, "mean"] <- 2 * .expr7/.expr9 .grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * .expr14)/.expr9^2) attr(.value, "gradient") <- .grad .value }
一応確認しておくと、
.expr2 <- sd^2 .expr7 <- x - mean .expr9 <- 2 * .expr2 .grad[, "mean"] <- 2 * .expr7/.expr9 # 2 * .expr7/.expr9 = # 2 * (x - mean) / 2 * sd^2 = # (x - mean) / sd^2
および、
.expr1 <- 2 * pi .expr2 <- sd^2 .expr3 <- .expr1 * .expr2 .expr7 <- x - mean .expr8 <- .expr7^2 .expr9 <- 2 * .expr2 .expr14 <- 2 * sd .grad[, "sd"] <- -(.expr1 * .expr14/.expr3/2 - .expr8 * (2 * .expr14)/.expr9^2) # -.expr1 * .expr14/.expr3/2 = # -2 * pi * 2 * sd / (2 * pi * sd^2) / 2 = # -1/sd # .expr8 * (2 * .expr14)/.expr9^2 = # (x - mean)^2 * (2 * 2 * sd) / (2 * sd^2)^2 = # (x - mean)^2 / sd^3
となります。
このぐらいなら手計算の方が早いですが、いろいろな場面で利用できるのではないでしょうか。
また、Hessian を自動で出してくれる deriv3 もあります。
注意
deriv のヘルプのタイトルは、
Symbolic and Algorithmic Derivatives of Simple Expressions
となっていて、以下の説明があります。
The internal code knows about the arithmetic operators +, -, *, / and ^, and the single-variable functions exp, log, sin, cos, tan, sinh, cosh, sqrt, pnorm, dnorm, asin, acos, atan, gamma, lgamma, digamma and trigamma, as well as psigamma for one or two arguments (but derivative only with respect to the first). (Note that only the standard normal distribution is considered.)
