パラメータ推定（最尤推定法）
今、手元にn個の観測値（標本）が得られているとする。
統計モデルとして表現するには、
が必要となる。
1番目の確率分布の形が分かったとして、2番目のパラメータを推定する方法の１つである「最尤推定法」に関するメモ。
ある1つの観測値が、他の観測値に影響を及ぼさない（確率変数は独立している、という）場合、それぞれの観測値の確率密度を掛けあわせたもの（同時確率）が尤度関数となる。
観測値つまり確率変数をx、未知パラメータをθとおくと、尤度関数は以下の形で表される。
\[
L(\theta) = \displaystyle \prod_{ i = 1 }^n f(x_i | \theta)
\]
尤度関数すなわち同時確率を最大にするパラメータθの値が「最尤推定値」になる。
小数の乗算はアンダーフローを引き起こすので、尤度関数の対数をとった「対数尤度関数」を最大化するパラメータθを求める。
具体的に見てみる。
確率分布の形が「正規分布」だとすると、推定したいパラメータθ は、
となる。
（μは平均、σは標準偏差を表す）
を対数尤度関数の式にあてはめて、式変形する。
導出過程は以下のサイトを参照。
統計学的画像再構成法である
これをRの関数として定義する。
第１項は定数なので、最大化を考えるにあたっては無視できるが、今回は確認用に残しておく。
対数尤度関数（正規分布）
# arg1 : vector(size 2), [1]mean, [2]variance # arg2 : vector(size N), observation log_likelihood_norm <- function(param, obs){ mu <- param[1]; sigma2 <- param[2]; term1 <- (-1) * length(obs) * 0.5 * log(2*pi) term2 <- (-1) * length(obs) * 0.5 * log(sigma2); term3 <- (-1) * 0.5 / sigma2 * sum((obs-mu)^2); return ( term1 + term2 + term3 ); }
検算。
平均「2」、分散「5」の正規分布に従う乱数を10個用意する。
> data <- rnorm(10, mean=2, sd=sqrt(5)) > data [1] -0.01731340 1.79861872 2.99598784 1.78515021 2.42412228 -1.61575147 [7] 1.74058535 3.84369888 1.93563703 0.03983286
Rのdnorm関数を使って求めた結果
value <- 0 for (i in 1:length(data)) { value <- value + log(dnorm(data[i], mean=2, sd=sqrt(5))) } value # -19.80805 value <- 1 for (i in 1:length(data)) { value <- value * dnorm(data[i], mean=2, sd=sqrt(5)) } log(value) # -19.80805
作成した関数での結果
log_likelihood_norm(c(2, 5), data) # -19.80805
一致していることを確認。
値が最大になるパラメータをRのoptim関数を使って探してみる。
最適解（サンプル数10）
optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.493105 2.313280 # $value # [1] -18.38531 # $counts # function gradient 73 NA # $convergence # [1] 0
出てきた最適値は、平均が「1.493105」、分散が「2.313280」
統計量（サンプル数10）
mean(data) # [1] 1.493057 # 不偏分散 var(data) # [1] 2.571645 # 標本分散 variance <- function(x) { var(x) * (length(x)-1) / length(x) } variance(data) # [1] 2.31448
サンプル数を100に増やしてみる。
最適解（サンプル数100）
data <- rnorm(100, mean=2, sd=sqrt(5)) optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.918190 4.851344 # $value # [1] -220.8375 # $counts # function gradient 59 NA # $convergence # [1] 0
出てきた最適値は、平均が「1.918190」、分散が「4.851344」
統計量（サンプル数100）
mean(data) # [1] 1.917897 var(data) # [1] 4.898471 variance(data) # [1] 4.849486
さらにサンプル数を増やしてみる。
最適解（サンプル数1000）
data <- rnorm(1000, mean=2, sd=sqrt(5)) optim(c(0.5, 1), fn=log_likelihood_norm, obs=data, control=list(fnscale=-1), hessian=FALSE) # $par # [1] 1.891001 4.752660 # $value # [1] -2198.32 # $counts # function gradient 63 NA # $convergence # [1] 0
統計量（サンプル数1000）
mean(data) # [1] 1.891252 var(data) # [1] 4.7577 variance(data) # [1] 4.752943
サンプル数が多くなるほど、標本分散と不偏分散の差が少なくなっている。
対数尤度関数が返す値の3次元散布図を書いてみる。
meanが「2の前後3、-1から5まで1刻みの7通り」、sigma2は「5の前後3、2から8まで1刻みの7通り」の49パターンとした。
一番高くなるところのパラメータが最尤推定値となる。
最尤推定値を求めるにあたって、Rのoptim関数では以下の最適化アルゴリズムが使える模様。
method="Nelder-Mead" # シンプレックス法、関数が微分できなくても可 method="BFGS" # 準ニュートン法 method="L-BFGS-B" # 準ニュートン法の変種(推定するパラメータのlower/upperを指定) method="CG" # 共役勾配法 method="SANN" # シミュレーテッド・アニーリング法、関数が微分できなくても可
methodパラメータの省略時デフォルトは「Nelder-Mead」（シンプレックス法）
optim function | R Documentation
The default method is an implementation of that of Nelder and Mead (1965), that uses only function values and is robust but relatively slow.
It will work reasonably well for non-differentiable functions.
（参考）多変量正規分布のパラメータ推定で、optim関数のmethodを変えて比較してみた結果
パラメータ推定（多変量正規分布） - ichou1のブログ
