I’m trying to practice Julia by coding an ordinal regression model. I have defined a log-likelihood function and it works as expected when I used Optim.jl
for optimization. To gain additional speed, I want to try using the autodiff features in Julia. However, they don’t seem to work with my log-likelihood function and I can’t figure out why. Any help is appreciated.
Log-Likelihood code
function neg_loglik(params::Vector, X::Matrix, Y_ord::Vector)
# Define the link function
F(x) = cdf(Normal(0,1), x)
# Number of categories
K = maximum(Y_ord)
# Define the subset of the params vector for the cutpoints and the betas
c = params[1:(K-1)]
β = params[K:end]
# Number of observations
n = length(Y_ord)
# Initialize a vector to store likelihood values
L = zeros(n)
# Return Inf negative log-likelihood if the cutpoints are out of order
if !issorted(c)
return(Inf)
end
for i in 1:n
# Current observation
y = Y_ord[i]
# Linear predictor
η = dot(X[i,:], β)
if y == 1
L[i] = F(c[1] - η)
elseif y == K
L[i] = 1 - F(c[K - 1] - η)
else
L[i] = F(c[y] - η) - F(c[y-1] - η)
end
end
ℓ = log.(L)
return(-sum(ℓ))
end
neg_loglik_c(params::Vector) = neg_loglik(params, X, Y_ord)
Trying autodiff
ForwardDiff.gradient(neg_loglik_c, vcat(sort!(rand(K)), zeros(p)))
Error message
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7})
Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
(::Type{T})(::T) where T<:Number at boot.jl:772
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
...
Stacktrace:
[1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7})
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}, i1::Int64)
@ Base ./array.jl:966
[3] neg_loglik(params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}, X::Matrix{Float64}, Y_ord::Vector{Int64})
@ Main ./code:36
[4] neg_loglik_c(params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}})
@ Main ./code:88
[5] chunk_mode_gradient(f::typeof(neg_loglik_c), x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:123
[6] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:21
[7] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(neg_loglik_c), Float64}, Float64, 7}}}) (repeats 2 times)
@ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/gradient.jl:17
[8] top-level scope
@ REPL[7]:1