Maximum Likelihood 2

Maximizing a Normal Likelihood

Công thức1

Hàm ‘log likelihood’ của một phân phối chuẩn có dạng như sau:

Hàm chứa 2

Dựa vào công thức trên ta viết hàm chứa hàm ‘log likelihood’ như sau:

make.NegLogLik <- function(data, fixed=c(FALSE,FALSE)) {
        params <- fixed
        function(p) {
                params[!fixed] <- p
                mu <- params[1]
                sigma <- params[2]
                a <- -0.5*length(data)*log(2*pi*sigma^2)
                b <- -0.5*sum((data-mu)^2) / (sigma^2)
                -(a + b)
        } 
}

Giống như bài trước, chúng ta tạo một mẫu có phân phối chuẩn với các tham số định trước và “giả vờ” không biết những giá trị này để đi tìm chúng.

set.seed(1); normals <- rnorm(100, 1, 2)

Chạy hàm chứa ở trên với đối số là mẫu vừa tạo.

nLL <- make.NegLogLik(normals)
nLL
## function(p) {
##                 params[!fixed] <- p
##                 mu <- params[1]
##                 sigma <- params[2]
##                 a <- -0.5*length(data)*log(2*pi*sigma^2)
##                 b <- -0.5*sum((data-mu)^2) / (sigma^2)
##                 -(a + b)
##         }
## <environment: 0xe35af40>
ls(environment(nLL))
## [1] "data"   "fixed"  "params"

Ước tính các tham số bằng hàm optim()

optim(c(mu = 0, sigma = 1), nLL)$par
##    mu sigma 
##   1.2   1.8

Cố định một tham số, tìm tham số còn lại

nLL <- make.NegLogLik(normals, c(FALSE, 2))
optimize(nLL, c(-1, 3))$minimum
## [1] 1.2
nLL <- make.NegLogLik(normals, c(1, FALSE))
optimize(nLL, c(1e-6, 10))$minimum
## [1] 1.8

Vẽ đồ thị biểu diễn

nLL <- make.NegLogLik(normals, c(1, FALSE))
x <- seq(1.7, 1.9, len = 100)
y <- sapply(x, nLL)
plot(x, exp(-(y - min(y))), type = "l")

plot of chunk unnamed-chunk-7

nLL <- make.NegLogLik(normals, c(FALSE, 2))
x <- seq(0.5, 1.5, len = 100)
y <- sapply(x, nLL)
plot(x, exp(-(y - min(y))), type = "l")

plot of chunk unnamed-chunk-8

Thư viện stats4

library(stats4)
set.seed(1); normals <- rnorm(100, 1, 2)
nll <- function(sigma, mu) {
                a <- -0.5*length(normals)*log(2*pi*sigma^2)
                b <- -0.5*sum((normals-mu)^2) / (sigma^2)
                return(-(a + b))
        } 
 
mleout <- mle(minuslogl=nll,start=list(sigma=2,mu=2))
summary(mleout)
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nll, start = list(sigma = 2, mu = 2))
## 
## Coefficients:
##       Estimate Std. Error
## sigma      1.8       0.13
## mu         1.2       0.18
## 
## -2 log L: 400
nll <- function(sigma, mu) {
                a <- -0.5*length(normals)*log(2*pi*sigma^2)
                b <- -0.5*sum((normals-mu)^2) / (sigma^2)
                return(-(a + b))
        } 
mleout <- mle(minuslogl=nll,start=list(sigma=2,mu=2),fixed=list(sigma=2))
summary(mleout)
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nll, start = list(sigma = 2, mu = 2), fixed = list(sigma = 2))
## 
## Coefficients:
##    Estimate Std. Error
## mu      1.2        0.2
## 
## -2 log L: 402
nll <- function(sigma, mu) {
                a <- -0.5*length(normals)*log(2*pi*sigma^2)
                b <- -0.5*sum((normals-mu)^2) / (sigma^2)
                return(-(a + b))
        } 
mleout <- mle(minuslogl=nll,start=list(sigma=2,mu=2),fixed=list(mu=1))
summary(mleout)
## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = nll, start = list(sigma = 2, mu = 2), fixed = list(mu = 1))
## 
## Coefficients:
##       Estimate Std. Error
## sigma      1.8       0.13
## 
## -2 log L: 401

  1. Thạm khảo tại Larry Wasserman (2003), “All of Statistics” ↩

  2. Xem thêm khóa R Programming, https://www.coursera.org/course/rprog ↩