Raftery copula #91
Raftery copula #91
Conversation
new file: test/test_plotly.jl
deleted: src/Plotly_functions.jl new file: test/RafteryTest.jl deleted: test/test_plotly.jl
modified: src/Copulas.jl modified: src/MiscellaneousCopulas/RafteryCopula.jl
There was a problem hiding this comment.
Thanks a lot for this one, I did not even know this copula !
Overall this looks good, I gave a few comments on the code, if you could go fix these that would be very nice. Also a test of the correctness of the pdf and cdf would be nice.
There is something that tickles me : the current state of your code, especially the R.d error, is showing that you did not run your own tests before sending this code. How are you working with this package ? I suspect that you are complicating your own life by not using the standard workflow, maybe we could discuss that together if you want ? Basically you should be doing the following :
1° In some julia environnement, type ] dev Copulas into the repl to developp the package.
2° Open ~/.julia/dev/Copulas.jl or the location of the package, setup this location to use my repo as upstream and your repo as origin
3° Make branches inside this repository directly and work there.
This will allow you to :
- Use directly your changes in any other julia environement by typing
]dev Copulas, especially usefull if you use Revise - Run the tests directly into your developping environement by
] testor using the vscode interface.
If you need some help with that i can help you
|
I work in vs code. Since I still don't have a good handle on Julia
packages, sometimes I need to make separate code and then join it together.
Abd… sometimes I delete the package on my computer and eat again to avoid.
I apologize if I'm making some mistakes.
El El jue, 30 nov 2023 a la(s) 5:56, Oskar Laverny ***@***.***>
escribió:
… ***@***.**** approved this pull request.
Thanks a lot for this one, I did not even know this copula !
Overall this looks good, I gave a few comments on the code, if you could
go fix these that would be very nice. Also a test of the correctness of the
pdf and cdf would be nice.
There is something that tickles me : the current state of your code,
especially the R.d error, is showing that you did not run your own tests
before sending this code. How are you working with this package ? I suspect
that you are complicating your own life by not using the standard workflow,
maybe we could discuss that together if you want ?
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> + RafteryCopula(d, θ)
+
+The Multivariate Raftery Copula of dimension d is arameterized by ``\\theta \\in [0,1]``
+
+```math
+C_{\\theta}(\\mathbf{u}) = u_{(1)} + \\frac{(1 - \\theta)(1 - d)}{1 - \\theta - d} \\left(\\prod_{j=1}^{d} u_j\\right)^{\\frac{1}{1-\\theta}} - \\sum_{i=2}^{d} \\frac{\\theta(1-\\theta)}{(1-\\theta-i)(2-\\theta-i)} \\left(\\prod_{j=1}^{i-1}u_{(j)}\\right)^{\\frac{1}{1-\\theta}}u_{(i)}^{\\frac{2-\\theta-i}{1-\\theta}}
+```
+
+
+where ``u_{(1)}, \\ldots , u_{(d)}`` denote the order statistics of ``u_1, \\ldots ,u_d``.
+
+More details about Multivariate Raftery Copula are found in :
+
+ Saali, T., M. Mesfioui, and A. Shabri, 2023: Multivariate Extension of Raftery Copula. Mathematics, 11, 414, https://doi.org/10.3390/math11020414.
+
+ Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.6.
this line should be :
***@***.***) Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.6.
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> +Constructor
+
+ RafteryCopula(d, θ)
+
+The Multivariate Raftery Copula of dimension d is arameterized by ``\\theta \\in [0,1]``
+
+```math
+C_{\\theta}(\\mathbf{u}) = u_{(1)} + \\frac{(1 - \\theta)(1 - d)}{1 - \\theta - d} \\left(\\prod_{j=1}^{d} u_j\\right)^{\\frac{1}{1-\\theta}} - \\sum_{i=2}^{d} \\frac{\\theta(1-\\theta)}{(1-\\theta-i)(2-\\theta-i)} \\left(\\prod_{j=1}^{i-1}u_{(j)}\\right)^{\\frac{1}{1-\\theta}}u_{(i)}^{\\frac{2-\\theta-i}{1-\\theta}}
+```
+
+
+where ``u_{(1)}, \\ldots , u_{(d)}`` denote the order statistics of ``u_1, \\ldots ,u_d``.
+
+More details about Multivariate Raftery Copula are found in :
+
+ Saali, T., M. Mesfioui, and A. Shabri, 2023: Multivariate Extension of Raftery Copula. Mathematics, 11, 414, https://doi.org/10.3390/math11020414.
Add the bibtex code to docs/src/assets/references.bib and use here the
following syntax :
***@***.***) Saali, T., M. Mesfioui, and A. Shabri, 2023: Multivariate Extension of Raftery Copula. Mathematics, 11, 414, https://doi.org/10.3390/math11020414.
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> +struct RafteryCopula{d, P} <: Copula{d}
+ θ::P # Copula parameter
+ function RafteryCopula(d,θ)
+ if (θ < 0) || (θ > 1)
+ throw(ArgumentError("Theta must be in [0,1]"))
+ elseif θ == 0
+ return IndependentCopula(d)
+ elseif θ == 1
+ return MCopula(d)
+ else
+ return new{d,typeof(θ)}(θ)
+ end
+ end
+end
+function _cdf(R::RafteryCopula, u::Vector{T}) where {T}
+ if length(u) != R.d
R.d does not exists here. You may access the dimension of the copula using
:
function _cdf(R::RafteryCopula{d,P}, u::Vector{T}) where {d,P,T}
or by using d = length(R), but R.d will not work.
Moreover, there is no need to check the dimension of u this is done
automatically now :)
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> + return MCopula(d)
+ else
+ return new{d,typeof(θ)}(θ)
+ end
+ end
+end
+function _cdf(R::RafteryCopula, u::Vector{T}) where {T}
+ if length(u) != R.d
+ throw(ArgumentError("Dimension mismatch"))
+ end
+
+ # Order the vector u
+ u_ordered = sort(u)
+
+ term1 = u_ordered[1]
+ term2 = (1 - R.θ) * (1 - R.d) / (1 - R.θ - R.d) * prod(u).^(1/(1 - R.θ))
R.d does not exist. check everywhere i will not mark it everytime
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> +
+function Distributions._rand!(rng::Distributions.AbstractRNG, R::RafteryCopula{d,P}, x::AbstractVector{T}) where {d,P,T <: Real}
+
+ dim = length(x)
+
+ # Step 1: Generate independent values u, u_1, ..., u_d from a uniform distribution [0, 1]
+ u = rand(rng, dim+1)
+
+ # Step 2: Generate j from a Bernoulli distribution with parameter θ
+ j = rand(Bernoulli(R.θ),1)
+ uj = u[1].^j
+ # Step 3: Calculate v_1, ..., v_d
+ for i in 2:dim+1
+ x[i] = u[i]^(1 - R.θ) * uj
+ end
+
x[1] was not set ? and x[d+1] was set, while x only has d values so this
will error.
------------------------------
In src/MiscellaneousCopulas/RafteryCopula.jl
<#91 (comment)>:
> + # Step 2: Generate j from a Bernoulli distribution with parameter θ
+ j = rand(Bernoulli(R.θ),1)
+ uj = u[1].^j
+ # Step 3: Calculate v_1, ..., v_d
+ for i in 2:dim+1
+ x[i] = u[i]^(1 - R.θ) * uj
+ end
+
+ return x
+end
+
+function ρ(R::RafteryCopula{d,R}) where {d, P}
+ term1 = (d+1)*(2^d-(2-R.θ)^d)-(2^d*R.θ*d)
+ term2 = (2-R.θ)^d*(2^d-d-1)
+ return term1/term2
+end
So Kendall tau would be good indeed, even if the formula is long, but most
importantly is the maximum likelyhood procedure that is developped in the
package : that would be nice to be able to fit this copula to some data !
------------------------------
In test/RafteryTest.jl
<#91 (comment)>:
> +
+ # Test PDF with some random values
+ u = rand(d)
+ pdf_value = pdf(F, u)
+ @test pdf_value >= 0
+ end
+end
+
***@***.*** "RafteryCopula Sampling" begin
+ using StableRNGs
+ rng = StableRNG(123)
+ n_samples = 100
+ F = RafteryCopula(3,0.5)
+ samples = rand(rng,F, n_samples)
+ @test size(samples) == (3, n_samples)
+end
Could you add you copula, for one or two parameter value, into the list in
test/margins_uniformity.jl ? this will automatically check that the
samples are uniform and that the cdf is OK with these samples. This is not
the best test, but it checks for coherence and usually catches a few bugs.
—
Reply to this email directly, view it on GitHub
<#91 (review)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/A7PTWNQJG2EPIZQXLUSYEJDYHBC3XAVCNFSM6AAAAABAAN4WZKVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTONJWHEYTENZWGY>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
|
No need to appologise, this is normal you have to learn :) Maybe the begining of this video could help you : https://www.youtube.com/watch?v=qM9NtiYlXck especially the |
|
Thanks for the tip!
I'm going to learn because I want to try to make my own package.
On the other hand, don't forget to tell me if you think it's better to put
Kendall's tau or not.
El El jue, 30 nov 2023 a la(s) 6:47, Oskar Laverny ***@***.***>
escribió:
… No need to appologise, this is normal you have to learn :)
Maybe this could help you : https://www.youtube.com/watch?v=qM9NtiYlXck
especially the Revise.jl part, if you dont use it yet
—
Reply to this email directly, view it on GitHub
<#91 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/A7PTWNQX3SSQYW5YARD7HRDYHBI37AVCNFSM6AAAAABAAN4WZKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMZTGQZDGMBTGE>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
|
Yes kendall tau would be nice ! |
|
Perfect, I'll work on it.
El El jue, 30 nov 2023 a la(s) 7:03, Oskar Laverny ***@***.***>
escribió:
… Yes kendall tau would be nice !
The fitting procedure from the article would be nice too !
—
Reply to this email directly, view it on GitHub
<#91 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/A7PTWNTTBQUGNIAZK54S3D3YHBKYDAVCNFSM6AAAAABAAN4WZKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMZTGQ2DONZYGU>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
|
The bibtex is next However, I don't know why the assert folder does not appear in my branch. I think I implemented the changes. Regarding the adjustment. It is simply a maximization of the likelihood function. However, when I use Distributions.fit it asks me for sufficient statistics and when I wanted to use optim to implement it independently it won't let me because of the dependencies. |
|
So I merged the main branch into yours, so that there is no more conflicts. Dont forget to pull the modifications or you will erase them ! For the fitting procedure, if it is just a non-linear maximization of the LLH, then let's not implement it. We should implement that somewhere else as a generic, I agree. |
Codecov ReportAttention:
Additional details and impacted files@@ Coverage Diff @@
## main #91 +/- ##
==========================================
- Coverage 80.70% 79.14% -1.56%
==========================================
Files 28 29 +1
Lines 679 729 +50
==========================================
+ Hits 548 577 +29
- Misses 131 152 +21 ☔ View full report in Codecov by Sentry. |
| term_numerator = 1 - d - R.θ * u_ordered[d]^((1 - R.θ - d) / (1 - R.θ)) | ||
| term_product = prod(u)^((R.θ) / (1 - R.θ)) | ||
|
|
||
| logpdf_value = log(term_numerator) - log(term_denominator) + log(term_product) |
There was a problem hiding this comment.
The pdf was not ran by the tests. Could you add a test to check if the pdf is, at least, not producing errors ? If possible, could you add another one that checks that the values produced are the right ones ? Maybe by computing one or two values by hand for specific parameters in dimension two ?
In particular, I am worried that term_denominator and term_numerator are negative, which will make the log produce an error. adding these tests is a very good way of knowing if it is ok or not :)
There was a problem hiding this comment.
I'm going to do it now
|
Perfect ! Let's merge that then, congratulation ! I will tag a new release to push it into the registry |
|
Great! Thank you for the suggestions
El El jue, 30 nov 2023 a la(s) 15:18, Oskar Laverny <
***@***.***> escribió:
… Merged #91 <#91> into main.
—
Reply to this email directly, view it on GitHub
<#91 (comment)>, or
unsubscribe
<https://github.com/notifications/unsubscribe-auth/A7PTWNV4KKWQADNVD7U3MXTYHDEW7AVCNFSM6AAAAABAAN4WZKVHI2DSMVQWIX3LMV45UABCJFZXG5LFIV3GK3TUJZXXI2LGNFRWC5DJN5XDWMJRGEYTENJXGA3DGOA>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
I added the multivariate Raftery copula and some tests. I did not implement Kendall's tau because it is a very long equation (its form is closed) however it is very long. If you think it is necessary, you could implement it. I leave you the main reference
[https://www.mdpi.com/2227-7390/11/2/414]
but you can also see details in
Nelsen 2006.