生物統計学

汎用最適化関数 optimize

optimize 関数はパラメーターを 1 つ持つ目的関数の最適化に用いられる。目的関数をおよび最大値・最小値の探索区間を代入することで簡単に利用できる。パラメーターを 2 つ以上持つ目的関数の最適化は optim 関数を利用する。

関数の最小値を求める

f(x) = (x - a)2 を最小にする x を求める例。ただし、a は定数として扱う。

x <- rnorm(10)
f <- function (x, a) {(x - a)^2}

opt <- optimize(f, interval = c(-5, 5), a = 2)
opt
## $minimum
## [1] 2
## $objective
## [1] 0

結果 opt オブジェクト中の minimum が、関数 f(x) を最小にする x である。また、 objective は関数 f(x) の最小値である。

探索区間を縮めると、以下のような結果を得ることができる。すなわち、探索区間内において、f(x) を最小にする x は 0.9999471 であり、このとき f(x = 0.9999471) = 1.000106 である。

optimize(f, interval = c(-5, 1), a = 2)
## $minimum
## [1] 0.9999471
## $objective
## [1] 1.000106

関数の最大値を求める

f(x) = -(x - a)2 を最大にする x を求める例。ただし、a は定数として扱う。この場合、optimize 関数の maximum 引数を TRUE にすればよい。 

x <- rnorm(10)
f <- function (x, a) {-(x - a)^2}

opt <- optimize(f, interval = c(-5, 5), a = 2, maximum = TRUE)
opt
## $maximum
## [1] 2
## $objective
## [1] 0

最小二乗法により平均を求める

正規分布に従う 10 個のデータから、それらの平均を最小二乗法により求める例を示す。10 個のデータをベクトル y に格納し、それらの平均を仮に x とおく。このような関数を求める

y <- rnorm(10, mean = 2)
mean(y)
## [1] 2.306972

f <- function(x, y) {sum((y - x)^2)}
optimize(f, interval = c(-10, 10), y = y)
## $minimum
## [1] 2.306972
## $objective
## [1] 4.347727

最尤法によりポアソン分布のパラメーターを求める

ポアソン分布はパラメーターを一つだけ持つ分布である。その確率関数は以下のように書ける。

\[ P(Y=y) = \frac{\lambda ^{y}}{y!}e^{-\lambda} \]

n 個の観測値がある場合、その対数尤度関数は以下のように書ける。

\[L(\lambda | \mathbf{y}) = \left(\sum_{i = 1}^{n}y_{i}\right)\log\lambda - \sum_{i=1}^{n}\log((y_{i})!) - n\lambda\]

ここで、ポアソン分布に従う乱数を 10 個生成し、それらのデータを利用してパラメーター λ を最尤法で求める。λ を最尤法で求めるには、対数尤度関数を最大にする λ を探せばよいので、以下のように実行できる。求めるパラメーター λ を変数 x と命名した。

R には階乗を求める関数がないため、階乗とガンマ関数の関係 n! = Gamma(n + 1) を利用した。

y <- rpois(10, lambda = 10)
mean(y)
## [1] 10.9

# 対数尤度関数
f <- function(x, y) {
  sum(y) * log(x) - sum(log(gamma(y + 1))) - length(y) * x
}
optimize(f, interval = c(0, 30), y = y, maximum = TRUE)
## $maximum
## [1] 10.9
## $objective
## [1] -24.76747

optimize で求めた結果は mean 関数と求めた結果と同じであることが確認できた。（ポアソン分布では E[Y] = Var[Y] = λ）