Условное максимальное правдоподобие для модели Раша

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

где — скрытая способность человека , а — трудность предмета .

Мы сымитируем данные на основе этой модели:

Random.seed!(123)
n = 1000
m = 5
theta = randn(n)
delta = randn(m)
r = zeros(n)
s = zeros(m)

for i in 1:n
  p = exp.(theta[i] .- delta) ./ (1.0 .+ exp.(theta[i] .- delta))
  for j in 1:m
    if rand() < p[j] ##correct
      r[i] += 1
      s[j] += 1
    end
  end
end
f = [sum(r.==j) for j in 1:m];

Поскольку число параметров растет с увеличением размера выборки, стандартное максимальное правдоподобие не даст нам состоятельных оценок. Вместо этого мы рассматриваем условное правдоподобие. Можно показать, что модель Раша является моделью экспоненциального семейства и что значение суммы является достаточной статистикой для . Если мы поставим условие на значение суммы, то сможем исключить . Действительно, с помощью алгебры мы можем показать следующее:

где — это элементарная симметричная функция порядка

где сумма вычисляется по всем возможным конфигурациям ответов, которые дают суммарное значение . Алгоритмы для вычисления и ее производных доступны в литературе (см., например, обзор от Бейкера (Baker) (1996) и более современный подход в работе Бискарри (Biscarri) (2018)).

function esf_sum!(S::AbstractArray{T,1}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  fill!(S,zero(T))
  S[1] = one(T)
  @inbounds for col in 1:n
    for r in 1:col
      row = col - r + 1
      S[row+1] = S[row+1] + x[col] * S[row]
    end
  end
end

function esf_ext!(S::AbstractArray{T,1}, H::AbstractArray{T,3}, x::AbstractArray{T,1}) where T <: Real
  n = length(x)
  esf_sum!(S, x)
  H[:,:,1] .= zero(T)
  H[:,:,2] .= one(T)

  @inbounds for i in 3:n+1
    for j in 1:n
      H[j,j,i] = S[i-1] - x[j] * H[j,j,i-1]
      for k in j+1:n
        H[k,j,i] = S[i-1] - ((x[j]+x[k])*H[k,j,i-1] + x[j]*x[k]*H[k,j,i-2])
        H[j,k,i] = H[k,j,i]
      end
    end
  end
end

esf_ext! (generic function with 1 method)

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

ϵ = ones(Float64, m)
β0 = zeros(Float64, m)
last_β = fill(NaN, m)
S = zeros(Float64, m+1)
H = zeros(Float64, m, m, m+1)

function calculate_common!(x, last_x)
  if x != last_x
    copyto!(last_x, x)
    ϵ .= exp.(-x)
    esf_ext!(S, H, ϵ)
  end
end
function neglogLC(β)
  calculate_common!(β, last_β)
  return -s'log.(ϵ) + f'log.(S[2:end])
end

neglogLC (generic function with 1 method)

Оценка параметров обычно выполняется относительно параметра без ограничений . Взяв производную по (и применив правило цепочки), получаем

где .

function g!(storage, β)
  calculate_common!(β, last_β)
  for j in 1:m
    storage[j] = s[j]
    for l in 1:m
      storage[j] -= ϵ[j] * f[l] * (H[j,j,l+1] / S[l+1])
    end
  end
end

g! (generic function with 1 method)

Аналогичным образом можно вычислить матрицу Гессиана

где .

function h!(storage, β)
  calculate_common!(β, last_β)
  for j in 1:m
    for k in 1:m
      storage[k,j] = 0.0
      for l in 1:m
        if j == k
          storage[j,j] += f[l] * (ϵ[j]*H[j,j,l+1] / S[l+1]) *
            (1 - ϵ[j]*H[j,j,l+1] / S[l+1])
        elseif k > j
          storage[k,j] += ϵ[j] * ϵ[k] * f[l] *
            ((H[k,j,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        else #k < j
          storage[k,j] += ϵ[j] * ϵ[k] * f[l] *
            ((H[j,k,l] / S[l+1]) - (H[j,j,l+1] * H[k,k,l+1]) / S[l+1] ^ 2)
        end
      end
    end
  end
end

h! (generic function with 1 method)

Оценки параметров элементов затем можно получить с помощью стандартных алгоритмов оптимизации (Ньютона-Рафсона или L-BFGS). Последняя проблема заключается в том, что модель не является идентифицируемой (умножение на константу и деление на ту же константу приводит к одинаковому правдоподобию). Поэтому при оценке параметров необходимо наложить какое-либо ограничение. Обычно либо \prod_{i=1}^m \epsilon_i = 1$ (что эквивалентно ).

con_c!(c, x) = (c[1] = sum(x); c)
function con_jacobian!(J, x)
  J[1,:] .= ones(length(x))
end
function con_h!(h, x, λ)
    for i in 1:size(h)[1]
        for j in 1:size(h)[2]
            h[i,j] += (i == j) ? λ[1] : 0.0
        end
    end
end
lx = Float64[]; ux = Float64[]
lc = [0.0]; uc = [0.0]
df = TwiceDifferentiable(neglogLC, g!, h!, β0)
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!, lx, ux, lc, uc)
res = optimize(df, dfc, β0, IPNewton())

 * Status: success

 * Candidate solution
    Final objective value:     1.157522e+03

 * Found with
    Algorithm:     Interior Point Newton

 * Convergence measures
    |x - x'|               = 1.25e-07 ≰ 0.0e+00
    |x - x'|/|x'|          = 4.76e-08 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.82e-13 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 5.89e-16 ≰ 0.0e+00
    |g(x)|                 = 3.46e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    23
    f(x) calls:    44
    ∇f(x) calls:   44

Сравните оценку с истинным положением дел.

delta_hat = res.minimizer
[delta delta_hat]

5×2 Matrix{Float64}:
 -0.915143  -0.387085
 -0.242781   0.242732
  1.22279    1.66615
 -3.05756   -2.62454
  0.667647   1.10274

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