Rによる数値微分 - 餡子付゛録゛

Rはderiv()関数で記号/微分ができるものの、和や積、つまりsum()やprod()があると微分できません。以下のようなコードを実行すると、『関数 'sum' は導関数の表中にありません』のようにエラーメッセージが出ます。

# ベクトルxを代入
x <- (1:10*23)%%7
# 変数bを2とする
b <- 2
# sum()を含む関数を設定
e <- expression(exp(b*sum(x)))
# 微分してみるとエラーが出る
deriv(e, "b")

解決策は調べた限りは三つあります。まずは微分してしまいましょう。

# 一階微分
sum(x)*exp(b*sum(x))
# 二階微分
(sum(x)^2)*exp(b*sum(x))

1.520372e+30と5.017226e+31が得られましたが、複雑な式になると便利ではありません。こういうときは数値微分を使います。
$f(x_0+h)$ を $x_0$ の周りでテイラー展開すると、以下のようになります。 $f(x_0+h)=f(x_0)+hf'(x_0)+\frac{h^2}{2}f''(x_0)+O(h^3)$
また $f(x_0-h)$ を $x_0$ の周りでテイラー展開すると、以下のようになります。
$f(x_0-h)=f(x_0)-hf'(x_0)+\frac{h^2}{2}f''(x_0)+O(h^3)$
ここで $f(x_0+h)-f(x_0-h)$ を整理すると、 $f'(x_0) = \frac{f(x_0 + h) - f(x_0 - h)}{2h} + O(h^2)$ となります（平野(2007)）。
これをRで表すと、

f <- function(b){
return(exp(b*sum(x)))
}
# 一階微分
df <- function(b, h=1e-4){
return(((f(b + h)-f(b - h))/(2*h)))
}
df(b)
# 二階微分
ddf <- function(b, h=1e-4){
return(((df(b + h)-df(b - h))/(2*h)))
}
ddf(b)

1.520372e+30と5.017226e+31が得られますが、手軽では無いかも知れません。
精度が欲しいときはnumDerivパッケージにgenD関数があるので、それを使いましょう。

library("numDeriv")
r <- genD(f, c(b))
r$D

一階微分と二階微分、両方を同時に求めてくれます。