パラメータ推定(最尤推定法)



今、手元にn個の観測値(標本)が得られているとする。

統計モデルとして表現するには、

  1. 観測値が従うであろう確率分布(確率密度関数)の形
  2. 確率密度関数のパラメータ

が必要となる。

1番目の確率分布の形が分かったとして、2番目のパラメータを推定する方法の1つである「最尤推定法」に関するメモ。

ある1つの観測値が、他の観測値に影響を及ぼさない(確率変数は独立している、という)場合、それぞれの観測値の確率密度を掛けあわせたもの(同時確率)が尤度関数となる。

観測値つまり確率変数をx、未知パラメータをθとおくと、尤度関数は以下の形で表される。
\[
L(\theta) = \displaystyle \prod_{ i = 1 }^n f(x_i | \theta)
\]

尤度関数すなわち同時確率を最大にするパラメータθの値が「最尤推定値」になる。
小数の乗算はアンダーフローを引き起こすので、尤度関数の対数をとった「対数尤度関数」を最大化するパラメータθを求める。

具体的に見てみる。
確率分布の形が「正規分布」だとすると、推定したいパラメータθ は、
 \theta = \left( \mu , \sigma^2 \right)
となる。
μは平均、σ標準偏差を表す)

正規分布確率密度関数


f(x | \mu, \sigma^2 )=\frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right)

を対数尤度関数の式にあてはめて、式変形する。


\begin{eqnarray}
\log L(\mu, \sigma^2 ) &=& \log \displaystyle \prod_{ i = 1 }^n f(x_i | \mu, \sigma^2 ) \\
 &= & - \frac{n}{2} \log 2\pi - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \displaystyle \sum_{ i = 1 }^{ n } (x_i - \mu)^2
\end{eqnarray}

導出過程は以下のサイトを参照。
統計学的画像再構成法である

これをRの関数として定義する。
第1項は定数なので、最大化を考えるにあたっては無視できるが、今回は確認用に残しておく。

対数尤度関数(正規分布
# 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パターンとした。

f:id:ichou1:20181017164527p:plain

一番高くなるところのパラメータが最尤推定値となる。

最尤推定値を求めるにあたって、Rのoptim関数では以下の最適化アルゴリズムが使える模様。
(methodパラメータの省略時デフォルトは「Nelder-Mead」(シンプレックス法))

method="Nelder-Mead"  # シンプレックス法、関数が微分できなくても可
method="BFGS"         # 準ニュートン法
method="L-BFGS-B"     # 準ニュートン法の変種(推定するパラメータのlower/upperを指定)
method="CG"           # 共役勾配法
method="SANN"         # シミュレーテッド・アニーリング法、関数が微分できなくても可

関数の最大・最小化 - RjpWiki