optimizationjuliapiecewiselog-likelihoodoptim

Optimization of piecewise functions in Julia


Extremely new to Julia, so please pardon any obvious oversights on my end

I am trying to estimate a piecewise likelihood function through optimization. I have the code functional in R, but have begun translating it to Julia in the hopes of faster estimation, for eventual bootstrapping

Here is the current block of code that I am trying (v and x are already as 1000x1 vectors elsewhere defined elsewhere):


function est(a,b)
    function pwll(v,x)
            if v>4 
                ILL=pdf(Poisson(exp(a+b*x)), v)
            elseif v==4
                ILL=pdf(Poisson(exp(a+b*x)), 4)+pdf(Poisson(exp(a+b*x)),3)+pdf(Poisson(exp(a+b*x)),2)
            else v==0
            ILL=pdf(Poisson(exp(a+b*x)), 1)+pdf(Poisson(exp(a+b*x)), 0)
            end
        return(ILL)
    end
    ILL=pwll.(v, x)
    function fixILL(x)
        if x==0
            x=0.00000000000000001
        else
            x=x
        end    
    end
    ILL=fixILL.(ILL)
    LILL=log10.(ILL)
    LL=-1*LILL
    return(sum(LL))
end


using Optim
params0=[1,1]
optimize(est, params0)

And the error message(s) I am getting are:

ERROR: InexactError: Int64(NaN)
Stacktrace:
 [1] Int64(x::Float64)
   @ Base ./float.jl:788
 [2] x_of_nans(x::Vector{Int64}, Tf::Type{Int64}) (repeats 2 times)
   @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/NLSolversBase.jl:60
 [3] NonDifferentiable(f::Function, x::Vector{Int64}, F::Int64; inplace::Bool)
   @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/nondifferentiable.jl:11
 [4] NonDifferentiable(f::Function, x::Vector{Int64}, F::Int64)
   @ NLSolversBase ~/.julia/packages/NLSolversBase/kavn7/src/objective_types/nondifferentiable.jl:10
 [5] promote_objtype(method::NelderMead{Optim.AffineSimplexer, Optim.AdaptiveParameters}, x::Vector{Int64}, autodiff::Symbol, inplace::Bool, args::Function)
   @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/optimize/interface.jl:63
 [6] optimize(f::Function, initial_x::Vector{Int64}; inplace::Bool, autodiff::Symbol, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/optimize/interface.jl:86
 [7] optimize(f::Function, initial_x::Vector{Int64})
   @ Optim ~/.julia/packages/Optim/tP8PJ/src/multivariate/optimize/interface.jl:83
 [8] top-level scope
   @ ~/Documents/Projects/ki_new/peicewise_ll.jl:120

I understand that it seems the error is coming from the function to be optimized being non-differentiable. A fairly direct translation works well in R, using the built in optim() function.

Can anyone provide any insight?

I have tried the above code displayed above, with multiple variations. The function to be optimized is functional, I am struggling with the optimization (the issues of which may stem from function being inefficiently written)


Solution

  • Here's an adapted version of your code which produces a solution:

    using Distributions, Optim
    
    function pwll(v, x, a, b)
        d = Poisson(exp(a+b*x))
        
        if v > 4 
            return pdf(d, v)
        elseif v == 4
            return pdf(d, 4) + pdf(d, 3) + pdf(d, 2)
        else
            return pdf(d, 1) + pdf(d, 0)
        end
    end
    
    fixILL(x) = iszero(x) ? 1e-17 : x
    
    est(a, b, v, x) = sum(-1 .* log10.(fixILL.(pwll.(v, x, a, b))))
    
    v = 4; x = 0.5 # Defining these here as they are not given in your post
    
    obj(input; v = v, x = x) = est(input[1], input[2], v, x)
    
    optimize(obj, [1.0, 1.0])
    

    I have no idea whether this is correct of course, check this against some sort of known result if you can.