Семинар по методу максимального правдоподобия

Наши данные подчиняются следующему закону:
\( x_i=i \), \( y_i=4.2 x_i+\varepsilon_i \), \( \varepsilon_i \sim N(0,9) \), \( i=1..100 \).

Создаём искуственно такие данные:

n.obs <- 100

x <- 1:n.obs
eps <- rnorm(n.obs, mean = 0, sd = 3)
y <- 4.2 * x + eps

А теперь сделаем вид, что мы забыли, чему равны истинные значения \( \beta \) и \( \sigma^2 \) и оценим модель
\( y_i=\beta x_i +\varepsilon_i \), где \( \varepsilon_i\sim N(0,\sigma^2) \):

log_lik <- function(beta, s2) {
    # это функция правдоподобия, она зависит от двух аргументов
    res <- -0.5 * log(s2) * n.obs - 0.5/s2 * sum((y - beta * x)^2)
    return(res)
}

neg_log_lik <- function(params) {
    # это функция правдоподобия с противоположным знаком она зависит от одного
    # векторного аргумента - так нужно для минимизации
    res <- (-1) * log_lik(params[1], params[2])
    return(res)
}

opt.res <- nlm(neg_log_lik, c(3, 2), hessian = T)
# с(3,2) - это стартовая точка. Лучшее её выбирать поближе к глобальному
# экстремуму. Если даже примерно не ясно, где находится глобальный
# экстремум, то попробуйте несколько разных стартовых значений, чтобы не
# попасться в локальный. Мы выбрали от фонаря.

opt.res$hessian  # гессиан в точке минимума

##           [,1]    [,2]
## [1,] 39245.175 -0.9560
## [2,]    -0.956  0.6724

opt.res$estimate  # точка минимума

## [1] 4.202 8.621

opt.res$minimum  # минимум функции

## [1] 157.7

var.hat <- solve(opt.res$hessian)  # оценка ковариационной матрицы


estimates <- opt.res$estimate  # вектор оценок неизвестных параметров
se <- sqrt(diag(var.hat))  # вектор ст. ошибок, корни из диагонали ковариационной матрицы


ci.left <- estimates - 1.96 * se
ci.right <- estimates + 1.96 * se

coef.table <- data.frame(estimates, se, ci.left, ci.right)
rownames(coef.table) <- c("beta", "sigma^2")
colnames(coef.table) <- c("Оценка", "Ст. ошибка", "Левая граница", 
    "Правая граница")

coef.table

##         Оценка Ст. ошибка Левая граница Правая граница
## beta     4.202   0.005048         4.192          4.212
## sigma^2  8.621   1.219522         6.231         11.012

Можно найти кучу примеров “MLE in R” в гугле

Замечания: