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

library("maxLik")

Задача 1. Взяв несколько пачек m&m-син исследователь Вениамин посчитал количества эмэндэмсин разных цветов.

x <- c(25, 35, 30, 32, 27, 31)

Вениамин хочет проверить гипотезу \(H_0\): все цвета встречаются равновероятно, с вероятностью \(1/6\).

Решение.

Задаём функцию правдподобия для данного случая: \[ L(p)=p_1^{x_1}\cdot p_2^{x_2}\cdot \ldots \cdot p_6^{x_6} \]

Логарифмируем, получаем \(\ell(p)=\sum_{i=1}^{6} x_i \ln(p_i)\).

ell <- function(p) {
  ans <- sum(x*log(p))
  return(ans)
}

В данной задаче можно руками без компа найти максимум правдоподобия :)

p_ml <- x/180
p_ml

## [1] 0.1388889 0.1944444 0.1666667 0.1777778 0.1500000 0.1722222

p_h0 <- rep(1/6, 6)
p_h0

## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

Значения правдоподобия в точке экстремума и в точке, соответствующей верной \(H_0\):

ell(p_ml)

## [1] -321.4425

ell(p_h0)

## [1] -322.5167

Считаем статистику отношения правдоподобия и находим критическое значение

LR <- 2*(ell(p_ml)-ell(p_h0))
LR

## [1] 2.148436

qchisq(0.95, df=5)

## [1] 11.0705

Вывод: \(H_0\) не отвергается

Если максимизировать в R, то надо явно выделить пять независимых вероятностей и создать отдельную функцию:

ell2 <- function(p5) {
  last_pro <- 1-sum(p5)
  all_p <- c(p5, last_pro)
  ans <- ell(all_p)
  return(ans)
}

Находим оценки с помощью R:

model <- maxLik(logLik = ell2, start = rep(1/6,5))
model$estimate

## [1] 0.1388889 0.1944444 0.1666667 0.1777778 0.1500000

Задача 2.

Величины \(X_1\), \(X_2\), , \(X_{100}\) независимы и нормальны \(N(\mu, \sigma^2)\). Для имеющейся выборки оказалось, что \(\sum X_i=200\), \(\sum X_i^2=1100\).

Найдите ML оценки \(\mu\), \(\sigma^2\).

Проверьте гипотезу \(H_0\): \(\mu=3\) и \(\sigma^2=10\)

Задаем функцию правдоподобия

ell <- function(pars) {
  mu <- pars[1]
  s2 <- pars[2]
  ans <- -50*log(s2)-
    (1100+100*mu^2-2*200*mu)/(2*s2)
  return(ans)
}

Для наглядности значения лог-правдоподобия в нескольких точках

ell(c(1,1))

## [1] -400

ell(c(2,7))

## [1] -147.2955

ell(c(100,100))

## [1] -5035.759

Максимизируем функцию правдоподобия в R:

model <- maxLik(logLik = ell, start=c(0,1))
pars_ml <- model$estimate
pars_ml

## [1] 2 7

Заводим вектор для \(H_0\):

pars_h0 <- c(3,10)
pars_h0

## [1]  3 10

LR-статистика и критические значения

LR <- 2*(ell(pars_ml)-ell(pars_h0))
LR

## [1] 15.66749

qchisq(0.95, df=2)

## [1] 5.991465

Вывод: \(H_0\) отвергается