引き続き，R に用意されている「数値計算を行う道具」を紹介する．
x の関数 f(x) について，f(x) = 0 を満たす解を求める方法をニュートン法という．
R では関数 uniroot() でニュートン法が行える．以下では 0 ≦ x ≦ 2 の範囲において，f(x) = x3 - 2 が 0 となる x の値を求めている．
f <- function (x) x^3 - 2 # (1) f(x) を定義
uniroot(f, c(0, 2)) # (2) 範囲をc(0, 2)で指定
# 結果の root に解が入っている：解は 1.259934 となっている
$root
[1] 1.259934
$f.root
[1] 6.088618e-05
$iter
[1] 5
$estim.prec
[1] 6.103516e-05
|
関数 polyroot() で多項式の解（根）を求めることが出来る．例えば，p(x) = 5 + 4x + x2 の根を求める場合は c(5, 4, 1) を与える．
polyroot(c(5, 4, 1)) # x^2 + 4x + 5 = 0 の解 [1] -2+1i -2-1i polyroot(c(1, 2, 1)) # x^2 + 2x + 1 = 0 の解 [1] -1-1.665335e-16i -1+1.665335e-16i # 重解の場合は同じ値が複数返される |
結果に小さな数値誤差が混じっているが，これを見えなくするには関数 round() を用いればよい．
round( polyroot(c(1, 2, 1)), digits=3 ) # x^2 + 2x + 1 = 0 の解 [1] -1+0i -1+0i |
数値計算では，数値微分は誤差が大きくなりやすいので細心の注意を払う必要があるのだが，R では関数を微分することが出来るので，あまり誤差のことを心配する必要はない．
f <- expression( a*x^4 ) # (1) 関数 f を定義（関数 expression() を用いること） D(f, "x") # (2) D(f, 微分する変数) で微分する a * (4 * x^3) |
高階微分する場合は以下のような関数 DD を定義すると便利である．
DD <- function(expr, name, order = 1) {
if(order < 1) stop("'order' must be >= 1")
if(order == 1) D(expr, name)
else DD(D(expr, name), name, order - 1)
}
DD(f, "x", 3) # D(数式, 微分する変数, 微分する回数)
a * (4 * (3 * (2 * x)))
|
結果の数式を使ってさらに計算を行う場合は関数 deriv() を用いればよい．
g <- deriv(~ exp(-x^2), "x", func=T) # (1) 関数 f = exp(-x^2) を微分
g(2) # (2) g(2) を実行すると f(2) と f'(2) が表示される
[1] 0.01831564
attr(,"gradient")
x
[1,] -0.07326256
a <- g(2) # 二次利用する（ f'(2) の値を使用する）場合は・・・
attr(a, "gradient") # こうすればよい
x
[1,] -0.07326256
str(a) # 興味のある方は class(a) などを実行されたし
atomic [1:1] 0.0183
- attr(*, "gradient")= num [1, 1] -0.0733
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr "x"
|
関数 deriv() を用いた計算例の一覧を示す．
|
コマンド |
機能 |
|
f <- deriv( ~ x^2, "x", func=T) |
x2 を x で微分したものを関数 f(x) に代入． |
|
g <- deriv( ~ x^2 *y , c("x","y"), func = TRUE) |
x2y を「 x で微分したもの」「 y で微分したもの」それぞれの関数を g(x,y) に代入． |
|
h <- deriv( ~ x^2 * y * z, c("x","y"), function(x, y, z=4){} ) |
x2yz を「 x で微分したもの」「 y で微分したもの」それぞれの関数を h(x,y,z=4) に代入． |
deriv() の引数に hessian=T を与えることで hessian を計算することも出来る．さらに，deriv(... , hessian=T) と同じ働きをする関数 deriv3() が用意されている．
数値積分は関数 integrate() で実行出来る．
f <- function(x) x^2 # (1) 積分する関数を定義 integrate(f, 0, 1) # (2) integrate(関数, 積分範囲の下限, 積分範囲の上限) 0.3333333 with absolute error < 3.7e-15 # 結果は 0.3333333 |
R に元々用意されている関数を積分する場合は，新たに関数を定義する必要はない．
integrate(sin, 0, pi) 2 with absolute error < 2.2e-14 |
さらに，無限大（Inf）を範囲に取ること（広義積分）も可能である．以下では正規分布の密度関数を -Inf （-∞） 〜 1.96 の範囲で積分している．
integrate(dnorm, -Inf, 1.96) # dnorm:正規分布の密度関数 0.9750021 with absolute error < 1.3e-06 |
パッケージ adapt 中の関数 adapt(次元, lower=積分の下限, upper=積分の上限,
functn=関数名) を用いれば多次元の数値積分を行うことが出来る．adapt() の引数
value に積分値が，relerr に相対誤差が表示される．また，関数の引数 eps で相対誤差の最小値を指定することも出来る．
library(adapt)
mynorm <- function(z) {
(1/sqrt(2*pi))^length(z) * exp(-0.5 * sum(z * z))
}
adapt(2, lo = c(-4,-4), up = c(4,4), functn = mynorm)
value relerr minpts lenwrk ifail
0.993709 0.0004751543 263 73 0
|
optim() はデフォルトでは最小化を行なうのだが，最大化をするには制御リスト control の fnscale に負の値を与えればよい (例 ： fnscale=-1) ．また，本来の目的関数の代わりに，値にマイナスを掛けた目的関数を用意する方法もある．
以下では，ベルヌーイ試行 10 回を行って成功が 3 回であった場合のパラメータ（ここでは確率） p を推定している．
LL <- function(p) {
choose(10,3)*p^(3)*(1-p)^(7) # 同時密度分布
}
optim(1, LL, control=list(fnscale=-1)) # 初期値 1
$par # 関数を最小にするパラメータ値
[1] 0.3
$value
[1] 0.2668279 # 最小にするパラメータ値を与えたときの関数の値
$counts
function gradient
58 NA
$convergence
[1] 0
$message
NULL
Warning message:
Nelder-Mead 法による一次元最適化は信頼できません: optimize を使用してください in:
|
警告メッセージが出ているが・・・，引数
method には全 5 種類のアルゴリズム（"Nelder-Mead", "BFGS",
"CG", "L-BFGS-B", "SANN"）を指定することが出来る．
関数の最小化を行なう場合は関数 nlm() を用いる．この関数は非線形回帰を行う際によく用いられる．以下では，ベルヌーイ試行 10 回を行って成功が 3 回であった場合のパラメータ（ここでは確率） p を推定している．
LL <- function(p) {
-choose(10,3)*p^(3)*(1-p)^(7) # 同時密度分布
}
nlm(LL, 0.5, fscale=-1) # 初期値 0.5
$minimum # 関数の最小値
[1] -0.2668279
$estimate # 関数が最小値をとるときのパラメータ値
[1] 0.3
$gradient
[1] 2.270406e-10
$code
[1] 1
$iterations
[1] 5
|
非線形回帰の詳しい解説は『工学のためのデータサイエンス入門』(間瀬
茂 他 著，数理工学社) を参照のこと．
パッケージ lpSolve の中の関数 lp() で数理計画を行うことが出来る．以下では条件（f.con）
x1 + 2x2 + 3x3 ≦ 9
3x1 + 2x2 + 2x3 ≦ 15
において，式（f.obj） x1 + 9x2 + x3 を最大化する例を挙げている．
library(lpSolve)
f.obj <- c(1, 9, 1)
f.con <- matrix (c(1, 2, 3, 3, 2, 2), nrow=2, byrow=TRUE)
f.dir <- c("<=", "<=")
f.rhs <- c(9, 15)
lp ("max", f.obj, f.con, f.dir, f.rhs)
Success: the objective function is 40.5
lp ("max", f.obj, f.con, f.dir, f.rhs)$solution
[1] 0.0 4.5 0.0
|
「変数が整数しかとらない」という条件を付ける場合は以下のようにする．
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3)
Success: the objective function is 37
lp ("max", f.obj, f.con, f.dir, f.rhs, int.vec=1:3)$solution
[1] 1 4 0
|
関数 fft() で高速フーリエ変換が行える．また，関数 convolve(x, y) で 2 つの数列の畳み込みが行える．
x <- 1:4 fft(x) [1] 10+0i -2+2i -2+0i -2-2i fft(fft(x), inverse = TRUE)/length(x) [1] 1+0i 2+0i 3+0i 4+0i |
関数
mvfff() も参照されたい．
ニュートン法