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

library(maxLik)

## Loading required package: miscTools

Реализуем ML своими руками через функцию 'optim'

Наши данные подчиняются следующему закону: \( 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) \):

Заметим, что параметр \( \beta \) произвольный, а параметр \( \sigma^2>0 \). Чтобы автоматически получать только положительные оценки параметра, часто используют простой трюк — перепараметризацию. То есть вводят новый параметр, равный логарифму неотрицательного параметра, \( a=\log(\sigma^2) \). Тогда \( \sigma^2 = \exp(a) \) будет автоматом положительным. Этот трюк позволяет избежать лишних ошибок при максимизации фукнции: мы гарантируем то, что алгоритм будет перебирать только положительные \( \sigma^2 \). 

m_log_lik <- function(params, y.in, x.in) {
    # это функция правдоподобия, она зависит от параметров и данных
    beta <- params[1]
    s2 <- exp(params[2])
    # этот трюк полезен, чтобы параметр s2 был заведомо неотрицательным

    res <- -0.5 * log(s2) * length(y.in) - 0.5/s2 * sum((y.in - beta * x.in)^2)
    # R умеет только минимизировать функции, поэтому возьмем функцию с минусом
    return(-res)
}

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

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

##            [,1]    [,2]
## [1,] 34843.5258 -0.8239
## [2,]    -0.8239 49.6466

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

## [1] 4.192 2.273

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

## [1] 163.3

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


estimates <- opt.res$par  # вектор оценок неизвестных параметров
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", "log(sigma^2)")
colnames(coef.table) <- c("Оценка", "Ст. ошибка", "Левая граница", 
    "Правая граница")

coef.table

##              Оценка Ст. ошибка Левая граница Правая граница
## beta          4.192   0.005357         4.182          4.203
## log(sigma^2)  2.273   0.141924         1.995          2.551

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

Реализуем ML с помощью пакета maxLik

log_lik <- function(params, y.in, x.in) {
    # это функция правдоподобия, она зависит от параметров и данных
    beta <- params[1]
    s2 <- exp(params[2])
    # этот трюк полезен, чтобы параметр s2 был заведомо неотрицательным

    res <- -0.5 * log(s2) * length(y.in) - 0.5/s2 * sum((y.in - beta * x.in)^2)
    return(res)
}

Оцениваем модель с помощью ML:

ml.results <- maxLik(log_lik, start = c(0, 0), y.in = y, x.in = x)

Выводим привычную табличку-отчёт:

summary(ml.results)

## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 13 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: -163.3 
## 2  free parameters
## Estimates:
##      Estimate Std. error t value Pr(> t)    
## [1,]  4.19211    0.00534     785  <2e-16 ***
## [2,]  2.26612    0.14139      16  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

Челябинские функции настолько суровы…

К сожалению, алгоритмы поиска экстремума еще далеки от совершенства. Есть функции, у которых человек легко найдёт экстремум устно, а популярные алгоритмы не смогут найти экстремум.

Например, \( f(x_1,x_2,x_3)=0.01\cdot (x_1-0.5)^2+ |x_1^2-x_2|+|x_1^2-x_3| \)

Если присмотреться, то легко заметить, что точка минимума \( x_1=0.5 \), \( x_2=0.25 \), \( x_3=0.25 \).

А что скажет компьютер?

bad.f <- function(x) {
    return(0.01 * (x[1] - 0.5)^2 + abs(x[1]^2 - x[2]) + abs(x[1]^2 - x[3]))
}
res <- optim(fn = bad.f, par = c(0, 0, 0))
res

## $par
## [1] 0.058622 0.003437 0.003432
## 
## $value
## [1] 0.001953
## 
## $counts
## function gradient 
##      501       NA 
## 
## $convergence
## [1] 1
## 
## $message
## NULL

Можно попробовать другие алгоритмы, встроенные в функцию optim.

Замечания: