丸めについて
R には実数を丸める以下の関数が用意されている．
|
関数 |
機能 |
|
round(x, digits = 0) |
IEEE 式で丸めを行なう（デフォルトは digits = 0 で小数以下を丸める）． |
|
ceiling(x) |
x 未満でない最小の整数を返す． |
|
floor(x) |
x 以上でない最大の整数を返す． |
|
signif(x, digits = 6) |
digits で指定された有効桁数(デフォルトは 6 桁)に丸める． |
|
trunc(x) |
x を「 0 へ向かって」整数化する． |
ニュートン法
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) に代入． |
数値積分
数値積分は関数 integrate() で実行出来る．
f <- function(x) x^2 # (1) 積分する関数を定義 integrate(f, 0, 1) # (2) integrate(関数, 積分範囲の下限, 積分範囲の上限) 0.3333333 with absolute error < 3.7e-15 # 結果は 0.3333333 |
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 |
関数の最大化・最尤推定法
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:
|
関数の最小化
関数の最小化を行なう場合は関数 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 |
Back to R
