#################################################
# B03: 參數估計                                 #
# 吳漢銘 國立政治大學統計學系                   #
# https://hmwu.idv.tw                           #
#################################################


# 7/19
y <- c(0.04304550, 0.50263474) #the observed sample
theta_hat <- length(y) / sum(y)
theta_hat

mlogL <- function(theta = 1) {
  n <- length(y)
  f <- -(n * log (theta) - theta * sum(y))
  f
}

library(stats4)
fit <- mle(mlogL)
summary(fit)



# 10/19
loglikefun <- function(x, par){
  mu <- par[1]
  sigma <- par[2]
  n <- length(x)
  # log of the normal likelihood
  # -n/2 * log(2*pi*s^2) + (-1/(2*s^2)) * sum((x-m)^2)
  loglikelihood <- - (n / 2)*(log(2 * pi * sigma^2)) + 
    (-1/(2 * sigma^2)) * sum((x - mu)^2)
  # return the negative to maximize rather than minimize
  - loglikelihood
}

set.seed(1123)
x <- rnorm(100)
x <- x/sd(x) * 8 # sd of 8
x <- x - mean(x) + 10 # mean of 10
cat("mean(x) =", mean(x), ", sd(x) =", sd(x))
optim(par = c(0.5, 0.5), fn = loglikefun, x = x)


# 12/19
alpha <- 0.05
xbar <- 21.2
sigma <- 8
n <- 100
v <- qnorm(1-alpha/2)*(sigma/sqrt(n))
c(xbar - v, xbar + v)