最尤推定の復習・その１

最尤推定法の復習

最尤推定法は、データに確率モデルを当てはめるときに必要になるパラメータを、データから推定する方法です。

私の場合、大学の教養課程で習ったものの在学中に使う機会はなく、就職してからも使う機会なく時間が流れ、40歳を目前としたころ急に必要になりました。習っていても使わぬまま20年過ぎてしまうと、覚えているはずもなく、Rではライブラリで処理できたことから、あいまいな知識のまま、今に至っている方法です。

私にとって最尤推定法を理解できない原因ははっきりしていて、常套句のように出てくる「観測されたデータをもっとも起こりやすい」という表現です。これがしっくりこない。白黒つかない印象がダメです。

とはいえ、これじゃあ、どうしようもないので、今の理解を記録しておこうと思います。

最尤推定法の手順

データ（観測値）を $X_1, X_2, \ldots, X_n$ とすると、尤度関数は $L(\theta) = f(X_1, X_2, \ldots, X_n | \theta)$ と書ける。この式が確率密度関数または確率関数になる。そして、この式の $\theta$ がパラメータで、これを推定するのが最尤推定法です。

また、最尤推定は $L(\theta) = f(X_1, X_2, \ldots, X_n | \theta)$ が成立するときの $\theta$ の最大値を求めるので、 $\hat{\theta} = \arg\max_{\theta} L(\theta)$ となり、この $\hat{\theta}$ を最尤推定値と呼ぶ。

そして、手順おおきく2つに分かれる。

1. 尤度関数を定義する データに当てはめたい数式が尤度関数になる。尤度関数にはパラメータが必要で、データが連続値なら確率密度関数、離散値なら確率関数のパラメータになる。
2. 尤度関数を最大化するパラメータの値を決める。 一般的に、最大化には最適化法を使うことが多く、最適化手法には様々な方法がある。

最尤推定の骨格

最尤推定法をRで実装するとき、最適化関数optimがある。最適化方法をオプションで指定できるので、コーディングするときに便利。

で、最尤推定法の手順で考えると、

尤度関数を定義する

パラメータとデータを引数にとるように尤度関数を定義する。

# 尤度関数
likelihood <- function(theta, data) {
  # 使いたい確率密度関数を書き込む
  # パラメータthetaとデータdataを引数にして、尤度を計算する
}

尤度関数を最大化するパラメータの値を決める

最適化方法にはいろいろあって、状況によって決まる。

パラメータの初期値を設定して、反復処理を繰り返して、最適値を探索する。

# 最尤推定
mle <- function(data) {
  # 尤度関数を最大化するパラメータの初期値を決める。
  par0 <- ...

  # 最適なパラメータを探索するための最適化手法を書き込む。
  # パラメータを最適化するためのパラメータ（ハイパーパラメータ）を
  # 定義することもある

  # 最大化した尤度と推定されたパラメータを返す。
  return(list(loglikelihood = ..., estimates = ...))
}

最適化方法はいくつかあるものの、optim関数を使えば、関数の定義が不要になるので、少し楽できる。

注意点としては、optim関数は尤度関数を最小化する関数なので、optim関数に渡す引数(fn)は最小化関数にする必要がある（対数尤度関数）。また、最適化手法やパラメータの初期値などをoptim関数の引数に指定するが、これは試行錯誤で決定する部分が多大にある。

# 最尤推定
mle <- function(data) {
  # optim関数は、デフォルトで尤度関数を最小化するので、
  # 対数尤度関数といった、尤度関数を最小化を考える。
  likelihood <- function(theta) {
    return(-likelihood(theta, data))
  }

  # 最適化関数optimを使ってパラメータを最適化する。
  # 尤度関数（fn）の他、パラメータの初期値（par）や最適化手法（method）も必要。
  result <- optim(par = ..., fn = likelihood, ...)

  # 最小化した負の対数尤度と推定されたパラメータを返す。
  return(list(loglikelihood = -result$value, estimates = result$par))
}