A mutable struct that does not get mutated ? or does it?

Specific Domains Statistics
1

Hi,

I’m trying to code an extension do Distributions.jl for a certain distribution (Some gamma convolutions), and implementing the Moschopoulos algorithm.

I have an issue: To compute the pdf, i need to cache certain values, that i compute on demand. When running pdf.(dist,1:10), the caching works correctly and each new call to pdf has values cached from the previous one.

On the other hand, when looking at the object dist after calling the function, cached values are gone.

Could someone explain to me what is going on ? Here is the code :


using Distributions
# Deal wich moshopoulos parameters :
mutable struct MoschopoulosParameters{T}
        θ₁::T
        C::T
        to_power::Array{T,1}
        γ::Array{T,1}
        ρ::T
        δ::Array{T,1}
end
function MoschopoulosParameters(α,θ)
    T = Base.promote_eltype(α,θ,[1.0])
    α = T.(α)
    θ = T.(θ)
    δ = T.([1])
    θ₁ = Base.minimum(θ)
    C = prod((θ ./ θ₁) .^ α)
    ρ = sum(α)
    to_power = (-θ₁ ./ θ) .+1
    γ = [sum(α .* to_power)] # gamma1
    return MoschopoulosParameters(θ₁,C,to_power,γ,ρ,δ)
end
eltype(P::MoschopoulosParameters) = typeof(P.θ₁)

struct UnivariateGammaConvolution{T<:Real} <: Distributions.ContinuousUnivariateDistribution where T
    α::Array{T}
    θ::Array{T}
    P::MoschopoulosParameters{T}
end

# Constructors :
UnivariateGammaConvolution(α::Real,θ::Real) = Distributions.Gamma(α,θ)
function UnivariateGammaConvolution(α::Array{T1,1},θ::Array{T2,1}) where {T1 <: Real, T2 <: Real}

    # Start by promoving alpha and theta to the same type :
    T = Base.promote_eltype(α,θ,[1.0]) # At least float.
    α = T.(α)
    θ = T.(θ)

    # Remove alphas and theta that are non-positives
    if any(α .== T(0))
        are_pos = α .* θ .!= T(0)
        α = α[are_pos]
        θ = θ[are_pos]
    end

    # regroup theta that are the sames :
    n = length(θ)
    for i in 1:(n-1)
        if i >= length(θ)
            break
        end
        for j in (i+1):length(θ)
            if θ[i] == θ[j]
                α[i] += α[j]
                deleteat!(α,j)
                deleteat!(θ,j)
                j = j-1
            end
        end
    end

    if length(α) == 1
        return Distributions.Gamma(α,θ)
    end
    P = MoschopoulosParameters(α,θ)
    tt = eltype(P)
    α = tt.(α)
    θ = tt.(θ)
    return UnivariateGammaConvolution(α,θ,MoschopoulosParameters(α,θ))
end

#### Support
Distributions.@distr_support UnivariateGammaConvolution 0.0 Inf

#### Conversions
eltype(d::UnivariateGammaConvolution) = typeof(d.α[1])
Base.convert(::Type{UnivariateGammaConvolution{T}}, d::UnivariateGammaConvolution{S}) where {T <: Real, S <: Real} = UnivariateGammaConvolution(T.(d.α), T.(d.θ))

#### parameters
shapes(d::UnivariateGammaConvolution) = d.α
scales(d::UnivariateGammaConvolution) = d.θ
rates(d::UnivariateGammaConvolution) = 1 / d.θ


#### Statistics
mean(d::UnivariateGammaConvolution) = sum(d.α .* d.θ)
var(d::UnivariateGammaConvolution) = sum(d.α .* d.θ .^ 2)


#### Sampling
function Base.rand(rng::Distributions.AbstractRNG,d::UnivariateGammaConvolution)
    sum(rand.(rng,Distributions.Gamma.(d.α,d.θ)))
end

#### Characteristic functions: pdf, cdf, mgf, cf

# Moshopoulos algorithm for pdf and cdf.
function MoschopoulosAlgo!(d::UnivariateGammaConvolution,x::Real, which, tolerence)

    @assert(which in ["pdf","cdf"], "which should be etiher pdf or cdf")
    print(length(d.P.γ),"at beginning\n")
    T = eltype(d)
    x = T(x)
    if x < T(0)
        return 0
    end
    k = 0
    out = T(0)
    while true
        if length(d.P.δ) < k+1
            # then compute the new deltas:
            n = length(d.P.δ)
            for k in n:(k+1)
                pushfirst!(d.P.γ, sum(d.α .* d.P.to_power .^ (k+1)))
                push!(d.P.δ,sum(d.P.γ[2:end] .* d.P.δ)/(k))
            end
        end
        dist = Distributions.Gamma(d.P.ρ + k,d.P.θ₁)
        if which == "pdf"
            step = d.P.δ[k+1] * Distributions.pdf(dist,x)
        elseif which == "cdf"
            step = d.P.δ[k+1] * Distributions.cdf(dist,x)
        end

        @assert(isfinite(step),"inf or nan append, the algorithm did not converge")
        out += step
        if isapprox(step,T(0),atol=tolerence,rtol=T(0))
            break
        end
        k = k+1
    end
    out *= d.P.C
    print(length(d.P.γ),"at end\n")
    return out
end

pdf(d::UnivariateGammaConvolution,x) = MoschopoulosAlgo!(d,x,"pdf", eps(eltype(d)))
cdf(d::UnivariateGammaConvolution, x) = MoschopoulosAlgo!(d,x,"cdf", eps(eltype(d)))
logpdf(d::UnivariateGammaConvolution, x) = log(pdf(d,x))
mgf(d::UnivariateGammaConvolution, t::Real) = prod((1 - t .* d.θ) .^ (-d.α))
cf(d::UnivariateGammaConvolution, t::Real) = prod((1 - (im * t) .* d.θ) .^ (-d.α))





# Let's now test that :
dist = UnivariateGammaConvolution([1,0.5],[4,2])
sample = zeros(Float64,100)
import Random
Random.rand!(dist,sample)
x = pdf.(dist,1:10)
print(length(dist.P.γ)) # Returns 1 !!

Edit : Looks like pdf(dist,x) caches the values correctly, while pdf.(dist,1:100) does not. I do not get why.

2

It probably has to do with the fact that you are broadcasting over dist.

By default, a dot call will try to iterate over all its arguments. Try wrapping dist in a Ref or a 1-element tuple: pdf.(Ref(dist),1:100) to prevent this from happening.

1 Like
3

It did not work, sadly. I dont get why a broadcasted call will not mutate it’s objects as a normal call, though.

Edit: The Ref() did not work, but the one ellement tuple did !

1 Like
4

I’m not sure exactly why, but since the broadcast call didn’t just entirely error out, I assume that your type UnivariateGammaConvolution must have a broadcasting rule defined implicitly for it by the Distributions.jl package because it is a subtype of Distributions.ContinuousUnivariateDistribution.

It might be worth opening an issue there to see if that is an intended behavior or not. I am also surprised that the Ref trick didn’t work, maybe mention that too.