library("maxLik")

Реализуем 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,] 37370.310266 -5.289489
## [2,]    -5.289489 49.834759
opt.res$par # точка минимума
## [1] 4.199285 2.203204
opt.res$value # минимум функции
## [1] 159.995
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.199285 0.005172969      4.189146       4.209424
## log(sigma^2) 2.203204 0.141656687      1.925557       2.480851

Можно найти кучу примеров “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: -159.9943 
## 2  free parameters
## Estimates:
##      Estimate Std. error t value Pr(> t)    
## [1,] 4.199143   0.005164  813.10  <2e-16 ***
## [2,] 2.199886   0.141390   15.56  <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.058621999 0.003436682 0.003431590
## 
## $value
## [1] 0.001953237
## 
## $counts
## function gradient 
##      501       NA 
## 
## $convergence
## [1] 1
## 
## $message
## NULL

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

Замечания: