StatsFuns.poislogccdf() and StatsFuns.logpdf() returning error message when used in a censored poisson regression

[sjwild]

I am trying a censored poisson regression using Turing, but unfortunately I am unable to do so. For context, I have some data where we asked respondents how many times they did an activity, but the options were top-coded. Naturally, this skews the responses.

When I try use Turing.@addlogprob! StatsFuns.poislogpdf() to try account for the censoring, I get the following error message:

ERROR: MethodError: no method matching poislogccdf(::ForwardDiff.Dual{ForwardDiff.Tag{Turing.Core.var"#f#3"{DynamicPPL.TypedVarInfo{NamedTuple{(:α, :β), Tuple{DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:α, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:α, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}, DynamicPPL.Metadata{Dict{AbstractPPL.VarName{:β, Tuple{}}, Int64}, Vector{Normal{Float64}}, Vector{AbstractPPL.VarName{:β, Tuple{}}}, Vector{Float64}, Vector{Set{DynamicPPL.Selector}}}}}, Float64}, DynamicPPL.Model{var"#51#52", (:y, :x, :censored, :N), (), (), Tuple{Vector{Int64}, Vector{Float64}, Vector{Int64}, Int64}, Tuple{}}, DynamicPPL.Sampler{NUTS{Turing.Core.ForwardDiffAD{40}, (), AdvancedHMC.DiagEuclideanMetric}}, DynamicPPL.DefaultContext}, Float64}, Float64, 2}, ::Int64)
Closest candidates are:
poislogccdf(::Union{Float64, Int64}, ::Union{Float64, Int64}) at /Users/stephenjwild/.julia/packages/StatsFuns/LvxaA/src/rmath.jl:86

Based on the message, I assume it is something to do with poislogccdf / poislogpdf within a Turing model. Using poislogccdf(exp(lambda[1]), y[1]) is fine if done outside the Turing model.

Here is some code for a minimum working example:

using Turing, Distributions, StatsFuns

# data
N_sims = 1000
x = rand(N_sims)
lambda = 0.3 .+ 0.5 .* x
y = rand.(Poisson.(exp.(lambda)))
censored = ifelse.(y .> 3, 1, 0)
y[y .> 3] .= 3

@model function poisson_censored(y, x, censored, N)

    # priors
    α ~ Normal(0, 2)
    β ~ Normal(0, 2)

    # model
    for n in 1:N
        if censored[n] == 1
            Turing.@addlogprob! StatsFuns.poislogccdf(exp.(α .+ β .* x[n]), y[n])
            #Turing.@addlogprob! logccdf.(Poisson(exp.(α .+ β .* x[n])), y[n])

        elseif censored[n] == 0
            Turing.@addlogprob! StatsFuns.poislogpdf(exp.(α .+ β .* x[n]), y[n])
            #Turing.@addlogprob! logpdf(Poisson(exp.(α .+ β .* x[n])), y[n])

        end

    end
    
end


mod_censored = poisson_censored(y, x, censored, N_sims)
chn_censored = sample(mod_censored, NUTS(), 2_000)


# works!
poislogccdf(exp(lambda[1]), y[1])

[devmotion]

The issue is not Turing specific. The problem is that StatsFuns.poslogccdf does not support dual numbers currently (it still uses the Rmath backend) (and it is also what logccdf(Poisson(), ...) ends up calling so it shouldn't matter what you use). JuliaStats/StatsFuns.jl#113 will improve this situation and replace it with a Julia implementation. AD should then work, possibly after defining some additional derivatives for specific functions in SpecialFunctions that are used by this implementation.

For the time being you can try to use the implementation in the PR instead of StatsFuns.poislogccdf.

[sjwild]

Awesome! Thanks for the response. Since this is already being addressed, I will close this issue.