頻度論とベイズ統計の主な考え方の違い 〜最尤法を理解する〜 [データサイエンス、統計モデル]
最近、ベイズ統計を教えたり自分で使ったりしていますが、最尤法の理解は避けて通れません。
そこで、最尤法の原理を理解するコードを作りました。
その前に、頻度論とベイズ統計の主な考え方の違いについて。
頻度論(ロナルド・フィッシャー)の世界
データが生成される確率分布のパラメータは点であると考える
(パラメータは定数で、データは確率分布)
ベイズ統計の世界
データが生成される確率分布を形成するパラメータも確率分布であると考える
(データは定数で、パラメータは確率分布)
ロナルド・フィッシャー(1890-1962)は、実験計画法、分散分析、最尤法といった現代の統計学の基礎を作った人です。
ところが、⼤のベイズ嫌いで、ベイズ統計は統計学ではないと言って駆逐しようとしました。w
続いて、最尤推定の理解について。
まず、パラメータは定数として、最初に神様が決めます。
その与えられたパラメータに従ってデータが観測されるのですが、確率分布として観測されるので、毎回同じ値が出てくるわけではないというのが基本思想となります。
例えば、平均0、分散1の正規分布からある値が出てきました。
その値が観測される確率は、
観測された値 0 → dnorm(0, 0, 1) = 0.3989423
観測された値 5 → dnorm(5, 0, 1) = 1.48672e-06
平均が0の正規分布なので、"0"が観測される確率は高く、0.398となっています。
"5"のような分布の端にある値が観測される確率は低く、1.48672e-06です。
今、平均と分散が分からない正規分布からデータが20個のデータを発生しました。
この時、平均と分散をどう推定したら良いか?
ちなみに、真の値ですが、
a <- 6
b <- 1
で設定しました。
そこで、最尤法の原理を理解するコードを作りました。
その前に、頻度論とベイズ統計の主な考え方の違いについて。
頻度論(ロナルド・フィッシャー)の世界
データが生成される確率分布のパラメータは点であると考える
(パラメータは定数で、データは確率分布)
ベイズ統計の世界
データが生成される確率分布を形成するパラメータも確率分布であると考える
(データは定数で、パラメータは確率分布)
ロナルド・フィッシャー(1890-1962)は、実験計画法、分散分析、最尤法といった現代の統計学の基礎を作った人です。
ところが、⼤のベイズ嫌いで、ベイズ統計は統計学ではないと言って駆逐しようとしました。w
続いて、最尤推定の理解について。
まず、パラメータは定数として、最初に神様が決めます。
その与えられたパラメータに従ってデータが観測されるのですが、確率分布として観測されるので、毎回同じ値が出てくるわけではないというのが基本思想となります。
例えば、平均0、分散1の正規分布からある値が出てきました。
その値が観測される確率は、
観測された値 0 → dnorm(0, 0, 1) = 0.3989423
観測された値 5 → dnorm(5, 0, 1) = 1.48672e-06
平均が0の正規分布なので、"0"が観測される確率は高く、0.398となっています。
"5"のような分布の端にある値が観測される確率は低く、1.48672e-06です。
今、平均と分散が分からない正規分布からデータが20個のデータを発生しました。
この時、平均と分散をどう推定したら良いか?
###################### 最尤推定の理解 ######################
# y = β0 + ε, ε~N(μ, σ2) から発生したモデルを推定する
# a <- ??? # 真の値
# b <- ??? # 真の値
# set.seed(1234)
# sample_dat <- round(rnorm(20, a, b), 3)
sample_dat <- c(6.415, 5.525, 6.066, 5.498, 5.174, 6.167, 5.104, 6.168, 6.355, 5.948,
5.804, 5.351, 4.890, 6.849, 6.022, 6.831, 4.756, 6.169, 6.673, 5.974)
sample_dat
# 自分で、a, bを推定する
a <- 0 # 仮に適当な数字を入れる
b <- 1 # 仮に適当な数字を入れる
# N(a, b)の時に、sample_dat(20レコード)が同時に観測される確率を計算する
# 尤度を計算
lhf <- dnorm(sample_dat, a, b)
lhf
prod(lhf) # 非常に小さい値
# 対数尤度を計算
sum(log(lhf))
# # 自分で対数尤度関数が最大となるa, bを見つける
a <- 0 # 仮に適当な数字を入れる
b <- 1 # 仮に適当な数字を入れる
# 対数尤度を計算
lhf <- dnorm(sample_dat, a, b)
sum(log(lhf))
ちなみに、真の値ですが、
a <- 6
b <- 1
で設定しました。