Extreme copulas

[Santymax98]

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]

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]

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 <#205 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7PTWNVWM4FLNIVV2RYE273ZH3WP7AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZTGU3DGNZXGY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[lrnv]

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 Distributions._logpdf 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

[Santymax98]

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 <#205 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7PTWNRER5LW466GNFLQBQTZH3XPFAVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNZTGU4DGNZVGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[lrnv]

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

Copulas.jl/src/Generator.jl

Lines 35 to 44 in cec1fab

	ϕ( G::Generator, t) = throw("This generator has not been defined correctly, the function `ϕ(G,t)` is not defined.")
	ϕ⁻¹( G::Generator, x) = Roots.find_zero(t -> ϕ(G,t) - x, (0.0, Inf))
	ϕ⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ(G,x), t)
	function ϕ⁽ᵏ⁾(G::Generator, k, t)
	X = TaylorSeries.Taylor1(eltype(t),k)
	taylor_expansion = ϕ(G,t+X)
	coef = TaylorSeries.getcoeff(taylor_expansion,k)
	der = coef * factorial(k)
	return der
	end

would help solve your differentiations issues ?

[Santymax98]

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

[lrnv]

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]

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]

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]

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]

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

[Santymax98]

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

[lrnv]

@Santymax98 I have reworked a bit your code, and I now have a different error :

julia> C = TEVCopula(4, [1.0 0.5; 0.5 1.0])
TEVCopula{2, 4, Matrix{Float64}}(
Σ: [1.0 0.5; 0.5 1.0]
)


julia> u=rand(2)
2-element Vector{Float64}:
 0.8992252448174132
 0.020655541515457676

julia> pdf(C,u)
ERROR: MethodError: no method matching _beta_inc(::ForwardDiff.Dual{…}, ::ForwardDiff.Dual{…}, ::ForwardDiff.Dual{…})
Stacktrace:
  [1] beta_inc(a::ForwardDiff.Dual{…}, b::ForwardDiff.Dual{…}, x::ForwardDiff.Dual{…})
    @ SpecialFunctions C:\Users\lrnv\.julia\packages\SpecialFunctions\QH8rV\src\beta_inc.jl:732
  [2] betaccdf
    @ C:\Users\lrnv\.julia\packages\StatsFuns\mrf0e\src\distrs\beta.jl:48 [inlined]
  [3] fdistccdf(ν1::ForwardDiff.Dual{…}, ν2::ForwardDiff.Dual{…}, x::ForwardDiff.Dual{…})
    @ StatsFuns C:\Users\lrnv\.julia\packages\StatsFuns\mrf0e\src\distrs\fdist.jl:18
  [4] tdistcdf(ν::ForwardDiff.Dual{…}, x::ForwardDiff.Dual{…})
    @ StatsFuns C:\Users\lrnv\.julia\packages\StatsFuns\mrf0e\src\distrs\tdist.jl:18
  [5] tdistcdf(ν::T, x::T) where T<:Real
    @ StatsFuns C:\Users\lrnv\.julia\packages\StatsFuns\mrf0e\src\distrs\tdist.jl:21 [inlined]
  [6] cdf
    @ C:\Users\lrnv\.julia\packages\Distributions\fgrZq\src\univariates.jl:647 [inlined]
  [7] ℓ(C::TEVCopula{2, 4, Matrix{Float64}}, u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{…}, ForwardDiff.Dual{…}, 2}})
    @ Copulas c:\Users\lrnv\.julia\dev\Copulas\src\ExtremeCopulas\tEVCopula.jl:49
  [8] #68
    @ c:\Users\lrnv\.julia\dev\Copulas\src\ExtremeValueCopula.jl:8 [inlined]
  [9] vector_mode_dual_eval!
    @ C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:24 [inlined]
 [10] vector_mode_gradient(f::Copulas.var"#68#70"{…}, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…}) 
    @ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\gradient.jl:89
 [11] gradient
    @ C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\gradient.jl:19 [inlined]
 [12] #99
    @ C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:16 [inlined]
 [13] vector_mode_dual_eval!
    @ C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:24 [inlined]
 [14] vector_mode_jacobian(f::ForwardDiff.var"#99#100"{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:125
 [15] jacobian(f::Function, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…}, ::Val{…})
    @ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:21
 [16] hessian(f::Function, x::Vector{…}, cfg::ForwardDiff.HessianConfig{…}, ::Val{…})
    @ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:17
 [17] hessian(f::Function, x::Vector{…}, cfg::ForwardDiff.HessianConfig{…})
    @ ForwardDiff C:\Users\lrnv\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:15
 [18] D_B_ℓ(C::TEVCopula{2, 4, Matrix{Float64}}, t::Vector{Float64}, B::Vector{Int64})
    @ Copulas c:\Users\lrnv\.julia\dev\Copulas\src\ExtremeValueCopula.jl:13
 [19] _logpdf(C::TEVCopula{2, 4, Matrix{Float64}}, u::Vector{Float64})
    @ Copulas c:\Users\lrnv\.julia\dev\Copulas\src\ExtremeValueCopula.jl:59
 [20] logpdf
    @ C:\Users\lrnv\.julia\packages\Distributions\fgrZq\src\common.jl:263 [inlined]
 [21] _pdf(d::TEVCopula{2, 4, Matrix{Float64}}, x::Vector{Float64})
    @ Distributions C:\Users\lrnv\.julia\packages\Distributions\fgrZq\src\common.jl:237
 [22] pdf(d::TEVCopula{2, 4, Matrix{Float64}}, x::Vector{Float64})
    @ Distributions C:\Users\lrnv\.julia\packages\Distributions\fgrZq\src\common.jl:222
 [23] top-level scope
    @ REPL[26]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia>

the first few lines shows that you are trying to differentiate the cdf of a multivariate T distribution, which is currently not implemented upstream. Maybe you could try to bypass this differentiation by writting things from the pdf and not the cdf ? I think that the pdf can be differentiated.

[lrnv]

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.

[Santymax98]

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 <#205 (review)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/A7PTWNQBX7UL6WSTWYNFWNLZIGP63AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMRYGYYTQMJUGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[lrnv]

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

[Santymax98]

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 > <#205 (review)>, > or unsubscribe > <https://github.com/notifications/unsubscribe-auth/A7PTWNQBX7UL6WSTWYNFWNLZIGP63AVCNFSM6AAAAABJODPZKOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDCMRYGYYTQMJUGM> > . > You are receiving this because you were mentioned.Message ID: > ***@***.***> >

[lrnv]

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