search

lrnv / Copulas.jl

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

Extreme copulas #205

Closed Santymax98 closed 2 months ago

Santymax98 commented 3 months ago

As you can see I have the comments in Spanish because I am still working on it. However, I gave some example so that you can see that the Galambos copula works well (The random number generation part is almost ready), however for the extreme value copula T I am having problems calculating the pdf. I made some changes so that it handles ForwarDiff.Dual but I'm still having problems

lrnv commented 3 months ago

Thanks @Santymax98. I will take a look and try to help you as soon as I have bandwidth. Do not hesitate to update this PR if/when/as-soon-as you make updates to it

Santymax98 commented 3 months ago

I will wait for your response. Until then, I will continue with other functions

El El lun, 17 jun. 2024 a la(s) 9:23, Oskar Laverny < @.***> escribió:

Thanks @Santymax98 https://github.com/Santymax98. I will take a look and try to help you as soon as I have bandwidth. Do not hesitate to update this PR if/when/as-soon-as you make updates to it

— Reply to this email directly, view it on GitHub https://github.com/lrnv/Copulas.jl/pull/205#issuecomment-2173563776, or unsubscribe https://github.com/notifications/unsubscribe-auth/A7PTWNVWM4FLNIVV2RYE273ZH3WP7AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZTGU3DGNZXGY . You are receiving this because you were mentioned.Message ID: @.***>

lrnv commented 3 months ago

A few remarks taking a glance at your code :

Santymax98 commented 3 months ago

Thanks for the tip!

Yes Oskar, clearly I will do the details as soon as I manage to solve the problems 🤣 also including some tests and verifications and obviously the documentation. I want everything to go in the same style that has been implemented so far.

El El lun, 17 jun. 2024 a la(s) 9:32, Oskar Laverny < @.***> escribió:

A few remarks taking a glance at your code :

  • Examples should go into the docs and test files and not at the end of the src/XXX/XXX.jl files
  • The _pdf and _cdf functions are internals that you have to implement but then in the examples you should use cdf() and pdf() from Distributions.jl : see in the src/Copula.jl file for the bindings between the two. So you should not export the _pdf,_cdf functions and rather use pdf,cdf in the tests.
  • using LinearAlgebra should go in the main src/Copulas.jl file.
  • Of course comments in english and correct documentation of the exported functions will be necessary later but i understand that you are not there yet :)
  • If you need to take recursives derivativs of a function take a look at how I implemented the derivatives of archimedean generators it may help you

— Reply to this email directly, view it on GitHub https://github.com/lrnv/Copulas.jl/pull/205#issuecomment-2173583753, or unsubscribe https://github.com/notifications/unsubscribe-auth/A7PTWNRER5LW466GNFLQBQTZH3XPFAVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZTGU4DGNZVGM . You are receiving this because you were mentioned.Message ID: @.***>

lrnv commented 3 months ago

Maybe implementing your derivatives like the generator derivates in src/Generator.jl:

https://github.com/lrnv/Copulas.jl/blob/cec1fab25ac7d3067b69d02b06b29a98b24fa91d/src/Generator.jl#L35-L44

would help solve your differentiations issues ?

Santymax98 commented 3 months ago

I think it changes a little because I need to calculate the mixed derivatives teorema

lrnv commented 3 months ago

Ok then. I think that your problem is that _pdf uses ForwardDiff in the ℓ function while ℓ(copula::TEVCopula, u::Vector) already uses frowardiff.

A good solution would be for ℓ(copula::TEVCopula, u::Vector) to not use forwarddiff: cant you compute these few derivatives by hand ? that is basically do the forwarddiff's job by hand in this function.

Santymax98 commented 3 months ago

I tried it and the derivative is very complicated, for the bivariate case it would have no problem (I even think it must be somewhere in the literature) but for d>2 calculating these derivatives requires a lot of work. For that reason I tried numerical methods. For Galambos it worked and for tEV it didn't. When I included ForwardDiff in ℓ(copula::TEVCopula, u::Vector) I had other types of errors such as that it is not possible to convert from dual to float64

lrnv commented 3 months ago

For Galambos it worked and for tEV it didn't

Yes, because ℓ(G::GalambosCopula, t::Vector) does not use ForwardDiff and does not over-type its internal variables.

The issue is your ℓ(copula::TEVCopula, u::Vector) function: it uses ForwardDiff and over-types its parameters and internal objects. Try untyping as much as you can in that function. Also convert_to_dual(Σ::Matrix{Float64}) explicitly types a lot of stuff, it might cause issues.

Santymax98 commented 3 months ago

If you change the line Σ = convert_to_dual(copula.Σ) to Σ = copula.Σ the types will be Float64 and will cause a conversion problem when calculating derivatives.

lrnv commented 3 months ago

Maybe you could use a multiply-by-one trick to promote to the right type ? like new_sigma = sigma * one(wanted_variable)

Santymax98 commented 3 months ago

I tried to rewrite as follows

    function ℓ(copula::TEVCopula, u::Vector)
    ν = copula.df
    Σ = copula.Σ
    d_len = length(u)
    l_u = zero(one(eltype(u)))  # Inicializar l_u con el tipo adecuado
    for j in 1:d_len
        R_j = construct_R_j(Σ, j)  # Construir la submatriz R_j
        println("R_j for j=$j: $R_j")  # Imprimir R_j para depuración

        terms = Vector{eltype(u)}()  # Asegurar que el vector terms tenga el tipo correcto
        for i in 1:d_len
            if i != j
                ρ_ij = Σ[i, j]
                term = (sqrt(ν + 1) * one(u[i]) / sqrt(one(u[i]) - ρ_ij^2)) * ((u[i] / u[j])^(-one(u[i])/ν) - ρ_ij * one(u[i]))
                push!(terms, term)
            end
        end
        println("terms for j=$j: $terms")  # Imprimir términos para depuración

        # Convertir temporalmente terms a Float64 para CDF
        terms_float64 = map(x -> ForwardDiff.value(x), terms)
        println("terms_float64 for j=$j: $terms_float64")

        if d_len == 2
            # Para el caso univariado, usar TDist con el parámetro adecuado
            T_dist = Distributions.TDist(ν + 1)
            par = terms_float64[1]
            println("TIPO DE PARAMETRO:", typeof(par))
            l_u += u[j] * Distributions.cdf(T_dist, par)
        else
            T_dist = Distributions.MvTDist(ν + 1, R_j)
            println("TIPOS DE PARAMETROS MULTIVARIADOS:", typeof.(terms_float64))
            l_u += u[j] * Distributions.cdf(T_dist, hcat(terms_float64...))
        end
    end
    return l_u
end

but different errors continue

Santymax98 commented 3 months ago

I also worked and found the error, in general my code was correct. The problem is that the T distribution, when you want to differentiate somewhere, uses the truncated beta function, which unfortunately is not compatible with dual-type variables. I implemented other copulas such as Husler Reiss and Marshal Olkin that have the same structure and did not cause problems.

El El mié, 19 jun. 2024 a la(s) 10:31, Oskar Laverny < @.***> escribió:

@.**** commented on this pull request.

So the error now has nothing to do with incorrect differentiation, just with the fact that the function you want to differentiate in-the-end calls a black box (probably some C code) that is not differentiable. I guess a correct approach would be to look for another, pure-julia, MvTDist implementation ? Maybe that exists somewhere.

— Reply to this email directly, view it on GitHub https://github.com/lrnv/Copulas.jl/pull/205#pullrequestreview-2128618143, or unsubscribe https://github.com/notifications/unsubscribe-auth/A7PTWNQBX7UL6WSTWYNFWNLZIGP63AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMRYGYYTQMJUGM . You are receiving this because you were mentioned.Message ID: @.***>

lrnv commented 3 months ago

So maybe this one shoudl be left aside. Looking forward for the other ones !

Santymax98 commented 3 months ago

I tried using SatatsFuns.tdist and the error was the same, I think a solution is to implement a TDist that supports dualForwardDiff

El El mié, 19 jun. 2024 a la(s) 10:32, Santiago Jimenez Ramos < @.***> escribió:

I also worked and found the error, in general my code was correct. The problem is that the T distribution, when you want to differentiate somewhere, uses the truncated beta function, which unfortunately is not compatible with dual-type variables. I implemented other copulas such as Husler Reiss and Marshal Olkin that have the same structure and did not cause problems.

El El mié, 19 jun. 2024 a la(s) 10:31, Oskar Laverny < @.***> escribió:

@.**** commented on this pull request.

So the error now has nothing to do with incorrect differentiation, just with the fact that the function you want to differentiate in-the-end calls a black box (probably some C code) that is not differentiable. I guess a correct approach would be to look for another, pure-julia, MvTDist implementation ? Maybe that exists somewhere.

— Reply to this email directly, view it on GitHub https://github.com/lrnv/Copulas.jl/pull/205#pullrequestreview-2128618143, or unsubscribe https://github.com/notifications/unsubscribe-auth/A7PTWNQBX7UL6WSTWYNFWNLZIGP63AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMRYGYYTQMJUGM . You are receiving this because you were mentioned.Message ID: @.***>

lrnv commented 3 months ago

Yes, that would be the right thing to do (but would take many efforts)