nlm はニュートン法を利用した最適化関数である。ニュートン法以外の最適化アルゴリズムも選択できる optim 関数の廉価版のような関数。

正規分布のパラメーター推定

正規分布のパラメーターとして平均と分散の 2 つある。optim 関数を利用して、最尤推定法によりパラメーターを推定する例を示す。まず、平均 10、分散 5 の正規分布に従う乱数を 100 個生成し、変数 y に格納する。

y <- rnorm(100, mean = 10, sd = 4)

次に正規分布のパラメーター推定時に用いる対数尤度関数を作成する。パラメーターを x とし、x[1] を平均、x[2] を分散とする。また、nlm 関数は最小化にしか対応していないため、対数尤度関数を最大化するために、対数尤度関数に -1 をかける。

loglikelihood <- function(x, y) {
  -1 * (- sum(0.5 * (y - x[1]) ^ 2 / x[2] + 0.5 * log(x[2])))
}

optim 関数を利用して対数尤度関数を最大にする x[1] と x[2] を求める。推定結果をみると、x[1] が 9.8、x[2] が 19.0 と乱数生成条件に近いことが確認できる。

nlm(loglikelihood, c(10, 2), y = y)
## $minimum
## [1] 197.133
## $estimate
## [1]  9.824181 18.966222
## $gradient
## [1] 7.782267e-07 1.078951e-06
## $code
## [1] 1
## $iterations
## [1] 9

負の二項分布のパラメーター推定

負の二項分布のパラメーターを最尤推定法により推定する例を示す。対数尤度関数の最大化は nlm 関数を利用する。尤度関数と確率分布関数が同じであることを利用して、ここでは尤度関数を自作せずに R で用意された負の二項分布の確率分布関数 dnbinom を利用する。nlm 関数は最小化を行う関数であるから、対数尤度関数に -1 をかける。

y <- rnbinom(100, mu = 1000, size = 200)

loglikelihood <- function(x, y){
  -1 * sum(log(dnbinom(y, mu = x[1], size = x[2])))
}

nlm 関数を利用してパラメーター推定を行う（対数尤度関数の最大化を行う）。

nlm(loglikelihood, c(900, 100), y = y)
## $minimum
## [1] 580.2502
## $estimate
## [1] 1001.3395  183.9273
## $gradient
## [1] 0 0
## $code
## [1] 1
## $iterations
## [1] 11