Оценка максимального правдоподобия: нормальная линейная модель

В следующем руководстве вы узнаете об оценке максимального правдоподобия в Julia для нормальной линейной модели.

Нормальная линейная модель (иногда называемая моделью OLS) является исполнительным компонентом регрессионного моделирования и используется во многих различных областях. В этом руководстве мы будем использовать смоделированные данные, чтобы продемонстрировать использование Julia для восстановления нужных параметров.

В первую очередь необходимо использовать пакет Optim, а также включить процедуру NLSolversBase:

using Optim, NLSolversBase
using LinearAlgebra: diag
using ForwardDiff

Добавьте Optim с помощью следующей команды в командной строке Julia: Pkg.add("Optim")

Сначала необходимо рассмотреть процесс создания данных или DGP. Следующий код выведет данные из нормальной линейной модели:

n = 40                              # Количество наблюдений
nvar = 2                            # Количество переменных
β = ones(nvar) * 3.0                # Истинные коэффициенты
x = [ 1.0   0.156651				# X-матрица объясняющих переменных плюс константа
 1.0  -1.34218
 1.0   0.238262
 1.0  -0.496572
 1.0   1.19352
 1.0   0.300229
 1.0   0.409127
 1.0  -0.88967
 1.0  -0.326052
 1.0  -1.74367
 1.0  -0.528113
 1.0   1.42612
 1.0  -1.08846
 1.0  -0.00972169
 1.0  -0.85543
 1.0   1.0301
 1.0   1.67595
 1.0  -0.152156
 1.0   0.26666
 1.0  -0.668618
 1.0  -0.36883
 1.0  -0.301392
 1.0   0.0667779
 1.0  -0.508801
 1.0  -0.352346
 1.0   0.288688
 1.0  -0.240577
 1.0  -0.997697
 1.0  -0.362264
 1.0   0.999308
 1.0  -1.28574
 1.0  -1.91253
 1.0   0.825156
 1.0  -0.136191
 1.0   1.79925
 1.0  -1.10438
 1.0   0.108481
 1.0   0.847916
 1.0   0.594971
 1.0   0.427909]

ε = [0.5539830489065279             # Ошибки
 -0.7981494315544392
  0.12994853889935182
  0.23315434715658184
 -0.1959788033050691
 -0.644463980478783
 -0.04055657880388486
 -0.33313251280917094
 -0.315407370840677
  0.32273952815870866
  0.56790436131181
  0.4189982390480762
 -0.0399623088796998
 -0.2900421677961449
 -0.21938513655749814
 -0.2521429229103657
  0.0006247891825243118
 -0.694977951759846
 -0.24108791530910414
  0.1919989647431539
  0.15632862280544485
 -0.16928298502504732
  0.08912288359190582
  0.0037707641031662006
 -0.016111044809837466
  0.01852191562589722
 -0.762541135294584
 -0.7204431774719634
 -0.04394527523005201
 -0.11956323865320413
 -0.6713329013627437
 -0.2339928433338628
 -0.6200532213195297
 -0.6192380993792371
  0.08834918731846135
 -0.5099307915921438
  0.41527207925609494
 -0.7130133329859893
 -0.531213372742777
 -0.09029672309221337]

y = x * β + ε;                      # Создать данные

В приведенном выше примере у нас есть 500 наблюдений, 2 объясняющие переменные плюс перехват, дисперсия ошибки равна 0,5, коэффициенты равны 3,0, и все они могут быть изменены пользователем. Поскольку мы знаем истинное значение этих параметров, мы должны получить эти значения при максимизации функции правдоподобия.

Следующим шагом будет определение функции Julia для функции правдоподобия. Следующая функция определяет функцию правдоподобия для нормальной линейной модели:

function Log_Likelihood(X, Y, β, log_σ)
    σ = exp(log_σ)
    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))
    llike = -llike
end

Log_Likelihood (generic function with 1 method)

Функция логарифмического правдоподобия принимает 4 входных значения: матрицу объясняющих переменных (X), зависимую переменную (Y), коэффициенты β и дисперсию ошибки. Обратите внимание, что во второй строке кода мы возводим в степень дисперсию ошибки, потому что она не может быть отрицательной, и хотим избежать этой ситуации при максимизации правдоподобия.

Следующим шагом является оптимизация функции. Сначала мы используем команду TwiceDifferentiable, чтобы в дальнейшем получить матрицу Гессиана, которая будет использоваться для формирования t-статистики:

func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),
                           ones(nvar+1); autodiff=:forward);

Приведенный выше оператор принимает 4 входных значения: матрицу x, зависимую переменную y, а также вектор β и дисперсию ошибки. vars[1:nvar] означает передачу вектора β, а vars[nvar + 1] — передачу дисперсии ошибки. Это можно рассматривать как вектор параметров, первые два из которых представляют β, а последний — дисперсию ошибки.

ones(nvar+1) являются начальными значениями параметров, а команда autodiff=:forward выполняет автоматическое дифференцирование в прямом режиме.

Собственно оптимизация функции правдоподобия выполняется с помощью следующей команды:

opt = optimize(func, ones(nvar+1))

 * Status: success

 * Candidate solution
    Final objective value:     1.722256e+01

 * Found with
    Algorithm:     Newton's Method

 * Convergence measures
    |x - x'|               = 2.79e-10 ≰ 0.0e+00
    |x - x'|/|x'|          = 9.12e-11 ≰ 0.0e+00
    |f(x) - f(x')|         = 0.00e+00 ≤ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00
    |g(x)|                 = 6.01e-14 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    9
    f(x) calls:    39
    ∇f(x) calls:   39
    ∇²f(x) calls:  9

Первым входным значением является функция, которую мы хотим оптимизировать, а вторым — начальные значения.

Через некоторое время вы увидите результаты работы процедуры оптимизации, причем оценки параметров будут очень близки к нашим смоделированным значениям.

Процедура оптимизации сохраняет несколько величин, и мы можем получить оценки максимального правдоподобия с помощью следующей команды:

parameters = Optim.minimizer(opt)

3-element Vector{Float64}:
  2.836642512088644
  3.053452125511052
 -0.9883745114256861

Имена полей для всех величин можно получить с помощью следующей команды: fieldnames(opt)

Чтобы получить правильную матрицу Гессиана, нам нужно «подставить» фактические значения параметров, которые максимизируют функцию правдоподобия, поскольку команда TwiceDifferentiable использует предпоследние значения для вычисления гессиана:

numerical_hessian = hessian!(func,parameters)

3×3 Matrix{Float64}:
 288.769        -21.7755        1.20224e-13
 -21.7755       223.84         -1.21827e-13
   1.20224e-13   -1.21827e-13  80.0

Найдем предполагаемое значение σ, а не логарифм σ, и его стандартную ошибку. Для этого воспользуемся дельта-методом: https://en.wikipedia.org/wiki/Delta_method.

эта функция экспоненциально увеличивает логарифм σ

function transform(parameters)
    parameters[end] = exp(parameters[end])
    parameters
end

transform (generic function with 1 method)

получение якобиана преобразования

J = ForwardDiff.jacobian(transform, parameters)'
parameters = transform(parameters)

3-element Vector{Float64}:
 2.836642512088644
 3.053452125511052
 0.3721811758462728

Теперь мы можем инвертировать матрицу Гессиана и использовать дельта-метод, чтобы получить матрицу дисперсий и ковариаций:

var_cov_matrix = J*inv(numerical_hessian)*J'

3×3 Matrix{Float64}:
  0.00348856   0.000339373  -1.75886e-18
  0.000339373  0.00450049    2.36094e-18
 -1.75886e-18  2.36094e-18   0.00173149

проверка корректности оценок параметров и t-статистики

t_stats = parameters./sqrt.(diag(var_cov_matrix))

3-element Vector{Float64}:
 48.02654897758057
 45.51568274694012
  8.94427190999916

просмотр результатов

println("parameter estimates:", parameters)
println("t-statsitics: ", t_stats)

parameter estimates:[2.836642512088644, 3.053452125511052, 0.3721811758462728]
t-statsitics: [48.02654897758057, 45.51568274694012, 8.94427190999916]

Далее можно ознакомиться с другими представляющими интерес статистическими данными, используя результаты работы программы оптимизации.

Простая программа

Ниже приводится версия программы без комментариев. Файл также доступен здесь: maxlikenlm.jl

using Optim, NLSolversBase
using LinearAlgebra: diag
using ForwardDiff

n = 40                              # Количество наблюдений
nvar = 2                            # Количество переменных
β = ones(nvar) * 3.0                # Истинные коэффициенты
x = [ 1.0   0.156651				# X-матрица объясняющих переменных плюс константа
 1.0  -1.34218
 1.0   0.238262
 1.0  -0.496572
 1.0   1.19352
 1.0   0.300229
 1.0   0.409127
 1.0  -0.88967
 1.0  -0.326052
 1.0  -1.74367
 1.0  -0.528113
 1.0   1.42612
 1.0  -1.08846
 1.0  -0.00972169
 1.0  -0.85543
 1.0   1.0301
 1.0   1.67595
 1.0  -0.152156
 1.0   0.26666
 1.0  -0.668618
 1.0  -0.36883
 1.0  -0.301392
 1.0   0.0667779
 1.0  -0.508801
 1.0  -0.352346
 1.0   0.288688
 1.0  -0.240577
 1.0  -0.997697
 1.0  -0.362264
 1.0   0.999308
 1.0  -1.28574
 1.0  -1.91253
 1.0   0.825156
 1.0  -0.136191
 1.0   1.79925
 1.0  -1.10438
 1.0   0.108481
 1.0   0.847916
 1.0   0.594971
 1.0   0.427909]

ε = [0.5539830489065279             # Ошибки
 -0.7981494315544392
  0.12994853889935182
  0.23315434715658184
 -0.1959788033050691
 -0.644463980478783
 -0.04055657880388486
 -0.33313251280917094
 -0.315407370840677
  0.32273952815870866
  0.56790436131181
  0.4189982390480762
 -0.0399623088796998
 -0.2900421677961449
 -0.21938513655749814
 -0.2521429229103657
  0.0006247891825243118
 -0.694977951759846
 -0.24108791530910414
  0.1919989647431539
  0.15632862280544485
 -0.16928298502504732
  0.08912288359190582
  0.0037707641031662006
 -0.016111044809837466
  0.01852191562589722
 -0.762541135294584
 -0.7204431774719634
 -0.04394527523005201
 -0.11956323865320413
 -0.6713329013627437
 -0.2339928433338628
 -0.6200532213195297
 -0.6192380993792371
  0.08834918731846135
 -0.5099307915921438
  0.41527207925609494
 -0.7130133329859893
 -0.531213372742777
 -0.09029672309221337]

y = x * β + ε;                      # Создать данные

function Log_Likelihood(X, Y, β, log_σ)
    σ = exp(log_σ)
    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))
    llike = -llike
end

func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),
                           ones(nvar+1); autodiff=:forward);

opt = optimize(func, ones(nvar+1))

parameters = Optim.minimizer(opt)

numerical_hessian = hessian!(func,parameters)

function transform(parameters)
    parameters[end] = exp(parameters[end])
    parameters
end

J = ForwardDiff.jacobian(transform, parameters)'
parameters = transform(parameters)

var_cov_matrix = J*inv(numerical_hessian)*J'

t_stats = parameters./sqrt.(diag(var_cov_matrix))

println("parameter estimates:", parameters)
println("t-statsitics: ", t_stats)

# Этот файл был создан с помощью Literate.jl, https://github.com/fredrikekre/Literate.jl.

Эта страница была создана с помощью Literate.jl.