ラボではうまくいったのに、本番はまったく。
ちなみに、基礎統計学はRをつかって教えているつもりです。
今日の授業に紹介したコードは最尤法についてです。
# 正規分布から尤度比を計算するための関数
llN = function(theta1, theta2, X) {
# This is the log-likelihood function for the normal distribution.
# X is the data (データのベクター)
# theta1 is the mean (i.e., average) (平均値)
# theta2 is the standard deviation (標準偏差)
# n is the length of the X vector (データの数)
# m is the length of the theta1 vector (推測したシータ１の値)
# pi is the circumference to diameter ratio (円周率)
n = length(X)
m = length(theta1)
if(length(theta1) == 1) {
# The log-likelihood function (尤度比関数)
L = -n/2*log(2*pi*theta2) + (-sum((X-theta1)^2) / (2*theta2))
} else {
L = rep(0, m)
for(i in 1:m) {
L[i] = -n/2*log(2*pi*theta2) + (-sum((X-theta1[i])^2) / (2*theta2))
}
}
return(L)
}
# 真の平均と標準偏差を指定
TRUE_MEAN = 2
TRUE_SD = 1
# シミュレーションデータを作成
DATA = rnorm(n = 50, mean = TRUE_MEAN, sd = TRUE_SD)
hist(DATA)
# DATAの平均
EST_MEAN = mean(DATA)
# 尤度比関数をつかってみよう。
# 平均を推定
my_guess = 3
llN(theta1 = my_guess , theta2 = TRUE_SD, X = DATA)
# 次は、推定平均の範囲を決めてから計算する。
guess = seq(from = -1, to = 3, length = 10)
L = llN(theta1 = guess, theta2 = 1, X = DATA)
# 尤度比を作図
plot(L ~ guess, type = "b", xlab = "guess", ylab = "log-likelihood")
abline(v = TRUE_MEAN, col = 2, lwd = 2)
abline(v = EST_MEAN, col = 3, lty = 2, lwd = 2)
legend("topleft", legend = c("log-likelihood","true mean","estimated mean"), col = c(1,2,3), lty = c(1,1,2))
