library("maxLik")

Задача. Величины \(X_1\), \(X_2\), , \(X_n\) независимы и имеют функцию плотности \(f(x)=\theta x^{\theta-1}\) при \(x\in [0;1]\). По 200 наблюдениям оказалось, что \(\sum \ln x_i=-220\).

Оцените параметр \(\theta\) методом максимального правдоподобия, проверьте \(H_0\): \(\theta=1\) с помощью трёх тестов.

Функция правдподобия равна \[ \ell(\theta)=200\ln \theta + (\theta-1)\cdot (-220) \]

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

ell <- function(theta) {
  ans <- 200*log(theta) + 
    (theta-1)*(-220)
  return(ans)
}

Пробуем функцию в нескольких точках

ell(0)

## [1] -Inf

ell(5)

## [1] -558.1124

Оцениваем модель:

report <- maxLik(logLik = ell, start = 1)
summary(report)

## --------------------------------------------
## Maximum Likelihood estimation
## Newton-Raphson maximisation, 4 iterations
## Return code 1: gradient close to zero
## Log-Likelihood: 0.937964 
## 1  free parameters
## Estimates:
##      Estimate Std. error t value Pr(> t)    
## [1,]  0.90909    0.06428   14.14  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## --------------------------------------------

В этой задаче экстремум можно найти руками без R

10/11

## [1] 0.9090909

Проверяем гипотезу \(H_0\): \(\theta=1\) с помощью LR-теста:

LR <- 2*(ell(10/11)-ell(1))
LR

## [1] 1.875928

Критическое значение равно:

qchisq(0.95, df=1)

## [1] 3.841459

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

Тест Вальда:

info <- summary(report)
se <- info$estimate[1,2]
(10/11-1)^2/se^2

## Std. error 
##   1.999972

Первая и вторая производные правдоподобия в точке \(\theta=1\) для LM-теста:

lp <- numericGradient(f=ell, t0 = 1)
lpp <- numericHessian(f=ell, t0 = 1)
lp

##      [,1]
## [1,]  -20

lpp

##           [,1]
## [1,] -199.9996

LM-статистика:

lp^2/(-lpp)

##          [,1]
## [1,] 2.000004