1

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)

Zachary Porreca
  • 13
  • 3
  • Can you make this an MWE? Your code doesn't run because `v` and `x` are undefined. – Nils Gudat Dec 14 '22 at 15:54
  • 2
    You may just need `params0=[1.0, 1.0]`. It looks like Optim is mutating this (or a copy of this) and trying to write a float into an array of integers is impossible. – mcabbott Dec 14 '22 at 16:02
  • Hi, I do not understand the line `else v==0` – after the else that is just a command – if v _is_ 0 this is just true but does nothing. Otherwise its false and also does nothing. – Ronny Dec 14 '22 at 16:03
  • Using @mcabbots hint – the `est` function (or cost for the optimiser) should be a function of one variable, not of 2. But I am not sure what you intend there, so that might be a next point to look for. – Ronny Dec 14 '22 at 16:06
  • Just saw your fix, you could define `add_eps(x) = [xi==0 ? eps(xi) : xi for xi in x]` and then call `add_eps(ILL)` (no dot) instead of your fix function. Looks maybe nicer – might be faster. – Ronny Dec 14 '22 at 16:09
  • Hey all, Thanks for these responses! @NilsGudat I'm not sure what an MWE is (very new to Julia), google keeps showing me links with the acronym-- from what I gather, its an executable of some sort? – Zachary Porreca Dec 14 '22 at 16:13
  • @mcabbott Still no luck after changing the arguments to params0=[1.0, 1.0] – Zachary Porreca Dec 14 '22 at 16:15
  • @Ronny Thanks for these suggestions. Would binding the a,b into a matrix (or Julia equivalent of that operation) allow for the single variable input? Where optimization output would be a 2x1 vector? The `else v==0` argument is necessary, meant to have the mass at both 0 and at 1 evaluated if the v input is equal to 0. Thanks for the "fix" fix. I appreciate that, and all of your help. Hopefully, I can get this running. – Zachary Porreca Dec 14 '22 at 16:20
  • MWE = Minimum Working Example. A piece of code that can be copy-pasted into a fresh Julia session and reproduces the error you are seeing. – Nils Gudat Dec 14 '22 at 16:42
  • I just do not know what the `else v==0` _means_, did you mean to write `elseif` ? though this is already a semantic problem. – Ronny Dec 14 '22 at 16:54
  • My two-input problem/question is answered in @nils-gudat s answer – to explain that: `Optim` always optimises _one_ variable and if I read the answer code correctly you have parameters `a` and `b` which are numbers and you want to optimise over? – Ronny Dec 14 '22 at 16:58

1 Answers1

3

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.

Nils Gudat
  • 13,222
  • 3
  • 39
  • 60