Merge pull request #37 from JoelTrent/dev

Update version to v0.2.0

Project.toml

@@ -1,7 +1,7 @@
 name = "EllipseSampling"
 uuid = "7d220c50-6298-48cd-9710-1d1a8ef13806"
 authors = ["JoelTrent <79883375+JoelTrent@users.noreply.github.com> and contributors"]
-version = "0.1.3"
+version = "0.2.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

src/generate_points.jl

@@ -382,7 +382,8 @@ function generate_N_clustered_points(num_points::Int, x_radius::T, y_radius::T,
 end
 
 """
-    generate_N_clustered_points(num_points::Int, Γ::Matrix{Float64}, θmle::Vector{Float64}, ind1::Int, ind2::Int; confidence_level::Float64=0.01, start_point_shift::Float64=rand(), sqrt_distortion::Float64=0.0)
+    generate_N_clustered_points(num_points::Int, Γ::Matrix{Float64}, θmle::Vector{Float64}, ind1::Int, ind2::Int; 
+        confidence_level::Float64=0.01, dof::Int=2, start_point_shift::Float64=rand(), sqrt_distortion::Float64=0.0)
 
 An alternate way to call [`generate_N_clustered_points(num_points::Int, e::Ellipse; start_point_shift::Float64=rand(), sqrt_distortion::Float64=0.)`](@ref), by supplying a square matrix Γ, the inverse of the Hessian of a log-likelihood function at its maximum likelihood estimate, indexes of the two variables of interest and the confidence level that represent a 2D ellipse approximation of the log-likelihood function.
 
@@ -395,13 +396,14 @@ An alternate way to call [`generate_N_clustered_points(num_points::Int, e::Ellip
 
 # Keyword Arguments
 - `confidence_level`: the confidence level ∈ [0.0,1.0] at which the ellipse approximation is constructed. Default is `0.01`.
+- `dof`: integer degrees of freedom used for calculation of the asymptotic confidence threshold defining the ellipse. Default is `2`.
 - `start_point_shift`: a number ∈ [0.0,1.0]. Default is `rand()` (defined on [0.0,1.0]), meaning that, by default, every time this function is called a different set of points will be generated.
 - `sqrt_distortion`: a number ∈ [0.0,1.0]. Default is `0.0`, meaning that, by default, this function will evenly space points on the the ellipse `e` with respect to the parameter `t`.
 """
 function generate_N_clustered_points(num_points::Int, Γ::Matrix{Float64}, θmle::Vector{Float64}, ind1::Int, ind2::Int; 
-    confidence_level::Float64=0.01, start_point_shift::Float64=rand(), sqrt_distortion::Float64=0.0)
+    confidence_level::Float64=0.01, dof::Int=2, start_point_shift::Float64=rand(), sqrt_distortion::Float64=0.0)
 
-    _, _, x_radius, y_radius, α = calculate_ellipse_parameters(Γ, ind1, ind2, confidence_level)
+    _, _, x_radius, y_radius, α = calculate_ellipse_parameters(Γ, ind1, ind2, confidence_level, dof)
     return generate_N_clustered_points(num_points, x_radius, y_radius, α, θmle[ind1], θmle[ind2],
                                     start_point_shift=start_point_shift, sqrt_distortion=sqrt_distortion)
 end

src/matrix_to_parameters.jl

@@ -1,14 +1,15 @@
 """
     calculate_ellipse_parameters(Γ::Matrix{Float64}, ind1::Int, ind2::Int,
-        confidence_level::Float64)
+        confidence_level::Float64, dof::Int)
 
-Given a square matrix Γ, the inverse of the Hessian of a log-likelihood function at its maximum likelihood estimate, indexes of the two variables of interest and the confidence level to construct a 2D ellipse approximation of the log-likelihood function, return the parameters of that ellipse.
+Given a square matrix Γ, the inverse of the Hessian of a log-likelihood function at its maximum likelihood estimate, indexes of the two variables of interest, the confidence level and degrees of freedom used to define ``\\ell_c``, which constructs a 2D ellipse approximation of the log-likelihood function, return the parameters of that ellipse.
 
 # Arguments
 - `Γ`: a square matrix which is the inverse of the Hessian of a log-likelihood function at its maximum likelihood estimate.
 - `ind1`: index of the first parameter of interest (corresponds to the row and column index of `Γ`)
 - `ind2`: index of the second parameter of interest (corresponds to the row and column index of `Γ`).
 - `confidence_level`: the confidence level ∈ [0.0,1.0] at which the ellipse approximation is constructed.
+- `dof`: integer degrees of freedom used for calculation of the asymptotic confidence threshold defining the ellipse. Default is `2`.
 
 # Details
 The parameters of interest are `a` and `b`, the radius of the major and minor axis respectively, `x_radius` and `y_radius`, the radius of the ellipse in the x and y axis respectively (i.e. the radius when the rotation `α` is zero) and `α`, an angle between 0 and π radians that the major axis of the ellipse has been rotated by from the positive x axis. `a` is equal to the maximum of `x_radius` and `y_radius`, while `b` is equal to the minimum of `x_radius` and `y_radius`.
@@ -35,23 +36,24 @@ where ``e_j`` and ``e_k`` are the ``j``th and ``k``th canonical vectors of ``\\m
 
 ## Obtaining Ellipse parameters
 
-By normalising our log-likelihood approximation equation for two parameters by our target confidence threshold of interest ``\\ell_c`` (at `confidence_level`) so that one side of the equation is equal to 1 we obtain the equation of an ellipse [friendlyelliptical2013](@cite):
+By normalising our log-likelihood approximation equation for two parameters by our target confidence threshold of interest ``\\ell_c`` (at `confidence_level`, with `dof` degrees of freedom, typically 2) so that one side of the equation is equal to 1 we obtain the equation of an ellipse [friendlyelliptical2013](@cite):
 ```math
 1 = -\\frac{1}{2\\ell_c} (\\psi-\\hat{\\psi})' ([e_j, e_k]' \\, \\Gamma(\\hat{\\theta}) \\, [e_j, e_k])^{-1} (\\psi-\\hat{\\psi}) = (\\psi-\\hat{\\psi})' \\mathcal{C} (\\psi-\\hat{\\psi}),
 ```
 
 The major and minor axis radii, `a` and `b` respectively, can then be evaluated by considering the inverse of the square roots of the eigenvalues of ``\\mathcal{C}`` (ordered from largest to smallest) [friendlyelliptical2013](@cite). To determine the rotation, `α` of the major axis of the ellipse from the positive ``x`` axis we calculate the inverse tangent of the division of the ``y`` and ``x`` components of the eigenvector corresponding to the largest eigenvalue [friendlyelliptical2013](@cite). 
 """
 function calculate_ellipse_parameters(Γ::Matrix{Float64}, ind1::Int, ind2::Int,
-    confidence_level::Float64)
+    confidence_level::Float64, dof::Int=2)
 
     (0.0 ≤ confidence_level && confidence_level ≤ 1.0) || throw(DomainError(confidence_level, "confidence_level must be between 0.0 and 1.0."))
+    (2 ≤ dof) || throw(DomainError(dof, "dof must be at least 2"))
     size(Γ)[1] == size(Γ)[2] || throw(DimensionMismatch("Γ must be a square matrix."))
     (0 < ind1 && ind1 ≤ size(Γ)[1]) || throw(BoundsError("ind1 must be a valid row index in Γ.")) 
     (0 < ind2 && ind2 ≤ size(Γ)[1]) || throw(BoundsError("ind2 must be a valid row index in Γ."))
 
     # normalise Hw so that the RHS of the ellipse equation == 1
-    Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level) * 0.5)
+    Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(dof), confidence_level) * 0.5)
     eigs = eigen(Hw)
     a_eig, b_eig = 1.0 ./ sqrt.(eigs.values)
 

test/runtests.jl

@@ -157,4 +157,44 @@ end
         @test isapprox_ellipsesampling(α, α_calc)
     end
 
+    @testset "CalculateEllipseParametersTest_HigherDof" begin
+        Γ = [7.7862e-6 -0.0506896 -0.0141446; -0.00506896 20.2146 6.61578; -0.0141446 6.61578 30.222]
+        ind1, ind2 = 2, 3
+        confidence_level = 0.01
+        dof=4
+        Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(dof), confidence_level)*0.5)
+        eigs = eigen(Hw)
+        a_eig, b_eig = sqrt.(1.0 ./ eigs.values)
+        eigvectors = eigs.vectors
+
+        a, b, x_radius, y_radius, α  = EllipseSampling.calculate_ellipse_parameters(Γ, ind1, ind2, confidence_level, dof)
+
+        @test isapprox_ellipsesampling(a_eig, a)
+        @test isapprox_ellipsesampling(b_eig, b)
+
+        α_eig = atan(eigvectors[2,1], eigvectors[1,1])
+        if α_eig < 0.0; α_eig += π end
+
+        @test isapprox_ellipsesampling(α_eig, α)
+
+        # For issue #30 when α → +/- 0.25 or +/- 1.25
+        a, b = 2.0, 1.0 
+        α = 0.25*π
+
+        Hw11 = (cos(α)^2 / a^2 + sin(α)^2 / b^2)
+        Hw22 = (sin(α)^2 / a^2 + cos(α)^2 / b^2)
+        Hw12 = cos(α)*sin(α)*(1/a^2 - 1/b^2)
+        Hw_norm = [Hw11 Hw12; Hw12 Hw22]
+
+        confidence_level=0.95
+        Hw = Hw_norm ./ (0.5 ./ (Distributions.quantile(Distributions.Chisq(dof), confidence_level)*0.5))
+        Γ = convert.(Float64, inv(BigFloat.(Hw, precision=64)))
+
+        a_calc, b_calc, _, _, α_calc = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level, dof)
+
+        @test isapprox_ellipsesampling(a, a_calc)
+        @test isapprox_ellipsesampling(b, b_calc)
+        @test isapprox_ellipsesampling(α, α_calc)
+    end
+
 end