lrnv / Copulas.jl

A fully `Distributions.jl`-compliant copula package
https://lrnv.github.io/Copulas.jl/
MIT License
83 stars 9 forks source link

[API choices] `fit_mle` for insteances, `_tau` #27

Open lrnv opened 1 year ago

lrnv commented 1 year ago
function Distributions.fit(D::SklarDist, x)
    # The first thing to do is to fit the marginals : 
    @assert length(D) == size(x, 1)
    # One could put fit but EM works for with fit_mle only
    m = Tuple(Distributions.fit_mle(D.m[i], x[i, :]) for i in axes(x, 1))
    u = pseudos(x)

    # τ⁻¹(FrankCopula, τ) works but not τ⁻¹(FrankCopula{3, Float64}, τ)
    # I found https://discourse.julialang.org/t/extract-type-name-only-from-parametric-type/14188/20
    # to just extract the typename.
    # Again defining fit(C::ArchimedeanCopula, u) = fit(::Type{CT},u) where {CT <: ArchimedeanCopula} would simplify the writting
    C = Distributions.fit(Base.typename(typeof(D.C)).wrapper, u)
    return SklarDist(C, m)
end

D̂mix = fit(D, simu)

Note that I had to define the weighted fit_mle(LogNormal, x, w) (only fit_mle(LogNormal, x) is in Distributions.jl). I believe this does the job

using StatsBase
function Distributions.fit_mle(::Type{<:LogNormal}, x::AbstractArray{T}, w::AbstractArray) where T<:Real
    lx = log.(x)
    μ, σ = mean_and_std(lx, weights(w))
    LogNormal(μ, σ)
end

(Note that StochasticEM() version only use unweighted version of fit_mle.)

Hence, IMO, it seems that for compatibility between EM.jl and Copula.jl, it is Copula.jl that needs to simply add Distributions.fit(D::SklarDist, x) as shown here. I believe it should not have any consequences on your existing functions.

Originally posted by @dmetivie in https://github.com/dmetivie/ExpectationMaximization.jl/issues/8#issuecomment-1481257615

lrnv commented 1 year ago

Maybe τ⁻¹(FrankCopula{3, Float64}, τ) should work and not only τ⁻¹(FrankCopula, τ) ? Dunno yet, have to check.