Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Bug] RafteryCopula fixes. #137

Merged
merged 5 commits into from
Feb 20, 2024
Merged

[Bug] RafteryCopula fixes. #137

merged 5 commits into from
Feb 20, 2024

Conversation

lrnv
Copy link
Owner

@lrnv lrnv commented Feb 12, 2024

Fixes #98

@Santymax98 it loks like the pdf is now fixed, but the cdf has values that still do not match what is in the tests files (see the @test_broken).

two month ago there #98 (comment) you proposed to use

function prueba(R::Vector{T}, u::Vector{T}) where T
    # Order the vector u
    u_ordered = sort(u)
    println("Sorted vector: ", u_ordered)
    
    term1 = u_ordered[1]
    println("Term 1: ", term1)
    
    term2 = ((1 - R[1]) * (1 - R[2])) / (1 - R[1] - R[2]) * prod(u)^(1/(1 - R[1]))
    println("Term 2: ", term2)
    
    term3 = 0.0
    for i in 2:length(R) # <<<<<<<<--------------- This should be 2:d and not 2:length(R) since length(R) is not the dimension. 
        prod_prev = prod(u_ordered[1:i-1])
        term3_part = (R[1] * (1 - R[1])) / ((1 - R[1] - i) * (2 - R[1] - i)) * prod_prev^(1/(1 - R[1])) * u_ordered[i]^((2 - R[1] - i) / (1 - R[1]))
        println("Term 3 (part $i): ", term3_part)
        term3 += term3_part
    end
    
    # Combine the terms to get the cumulative distribution function
    cdf_value = term1 + term2 - term3
    println("Final CDF value: ", cdf_value)
    
    return cdf_value
end

prueba([0.5,3], [0.1,0.2,0.3])

to check that the value's returned by the code are indeed correct. But your code has a mistake.

It looks to me like the current implementation is mathcing what is in the article, so maybe the issue is that the tests values are now wrong and should be updated ?

@lrnv lrnv changed the title Fixing the pdf and clearing out the code. [Bug] RafteryCopula fixes. Feb 12, 2024
@Santymax98
Copy link
Contributor

Santymax98 commented Feb 12, 2024 via email

@lrnv
Copy link
Owner Author

lrnv commented Feb 12, 2024

I remember that I did the tests with paper and pencil according to the formula in the document.

This is indeed a very good way of insuring that the implementation is OK. Maybe you made a mistake doing this, could you simply re-do it for the trivariate cdf tests ? the only ones that are still broken.

@Santymax98
Copy link
Contributor

function prueba_CDF(R::Vector{T}, u::Vector{T}) where T
    # Order the vector u
    θ = R[1]
    println("param:", θ)
    d = round(Int, R[2])
    println("dimension:", d)
    u_ordered = sort(u)
    println("Sorted vector: ", u_ordered)
    
    term1 = u_ordered[1]
    println("Term 1: ", term1)
    
    term2 = ((1 - θ) * (1 -d)) / (1 - θ - d) * prod(u)^(1/(1 - θ))
    println("Term 2: ", term2)
    
    term3 = 0.0
    for i in 2:d # <<<<<<<<--------------- This should be 2:d and not 2:length(R) since length(R) is not the dimension. 
        prod_prev = prod(u_ordered[1:i-1])
        term3_part = ((θ * (1 - θ)) / ((1 - θ - i) * (2 - θ - i))) * prod_prev^(1/(1 - θ)) * u_ordered[i]^((2 - θ - i) / (1 - θ))
        println("Term 3 (part $i): ", term3_part)
        term3 += term3_part
    end
    
    # Combine the terms to get the cumulative distribution function
    cdf_value = term1 + term2 - term3
    println("Final CDF value: ", cdf_value)
    
    return cdf_value
end

prueba_CDF([0.5,3], [0.1,0.2,0.3])

should be 0.08236... I calculated this value manually and the test code works fine. However, when I combine it with Copulas.jl it gives me the error that we have had so far.

@Santymax98
Copy link
Contributor

Santymax98 commented Feb 20, 2024

function prueba_PDF(R::Vector{T}, u::Vector{T}) where T
    # Order the vector u
    θ = R[1]
    d = round(Int, R[2])
    u_ordered = sort(u)
    println("Sorted vector: ", u_ordered)
    
    term1 = (1/(((1-θ)^(d-1))*(1-θ-d)))
    println("Term 1: ", term1)
    
    term2 = (1-d-θ*(u_ordered[d])^((1-θ-d)/(1-θ)))
    println("Term 2: ", term2)
    
    term3 = (prod(u))^((θ)/(1-θ))
    println(term3)
    # Combine the terms to get the density distribution function
    pdf_value = term1*term2*term3
    println("Final PDF value: ", pdf_value)
    
    return pdf_value
end
prueba_PDF([0.5,3], [0.1,0.2,0.3])

Regarding this test, the value in this case should be 1.99450... and on my machine the results are coming out correct.
Do you think that regarding density there is any problem with calculating the logarithm? Is it simpler to do it with multiplication?

@lrnv lrnv merged commit d8fed1e into main Feb 20, 2024
4 of 6 checks passed
@lrnv lrnv deleted the lrnv/issue98 branch February 20, 2024 09:17
@lrnv
Copy link
Owner Author

lrnv commented Feb 20, 2024

@Santymax98 I cleared up the code and fixed the previous value. I included your versions as fallbak testing. Thanks a lot this is now Ok for me !

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[Bug] RafteryCopula broken tests
2 participants