From 96be7ffd3e185ff991297e083cb47da45ba09283 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Mon, 22 Jan 2024 13:30:09 +1300
Subject: [PATCH 1/7] Add specification of the degrees of freedom, `dof`, used
 to define the boundary of the ellipse approximation of a log-likelihood
 function

---
 src/generate_points.jl      |  5 +++--
 src/matrix_to_parameters.jl | 12 +++++++-----
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/generate_points.jl b/src/generate_points.jl
index 9b7e9ef..ed33a2a 100644
--- a/src/generate_points.jl
+++ b/src/generate_points.jl
@@ -395,13 +395,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
\ No newline at end of file
diff --git a/src/matrix_to_parameters.jl b/src/matrix_to_parameters.jl
index 8f0f5c7..1568b7d 100644
--- a/src/matrix_to_parameters.jl
+++ b/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,7 +36,7 @@ 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}),
 ```
@@ -43,15 +44,16 @@ By normalising our log-likelihood approximation equation for two parameters by o
 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)
 

From bb2aa090a9977c9aa66a72cb98026633a343ab20 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Mon, 22 Jan 2024 13:33:11 +1300
Subject: [PATCH 2/7] Add `dof` kwarg to `generate_N_clustered_points`
 documentation definition

---
 src/generate_points.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/generate_points.jl b/src/generate_points.jl
index ed33a2a..017d37c 100644
--- a/src/generate_points.jl
+++ b/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.
 

From ab9076f84637ded9750a0808c3b4fa57141a9725 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Mon, 22 Jan 2024 13:53:05 +1300
Subject: [PATCH 3/7] Update version to v0.1.4

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 34afa31..c9f3c07 100644
--- a/Project.toml
+++ b/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.1.4"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

From 2e8ac3eefca27690c39fc29c17b7a5d31ad9c4ff Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Mon, 22 Jan 2024 14:05:27 +1300
Subject: [PATCH 4/7] Update version to v0.2.0 instead of v0.1.4, as api
 updated

---
 Project.toml     | 2 +-
 docs/make.jl     | 1 +
 test/runtests.jl | 1 +
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index c9f3c07..9ab64e9 100644
--- a/Project.toml
+++ b/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.4"
+version = "0.2.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/docs/make.jl b/docs/make.jl
index 10a5400..81305e5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,3 +1,4 @@
+using Revise
 using EllipseSampling
 using Documenter, DocumenterCitations
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 2bdb38b..a0f6ab9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,3 +1,4 @@
+using Revise
 using EllipseSampling
 using Test
 using LinearAlgebra

From 6028d5162daba1c4bcd66a7bbc4484957dcf43b2 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Tue, 23 Jan 2024 10:32:24 +1300
Subject: [PATCH 5/7] Revert "Update version to v0.2.0 instead of v0.1.4, as
 api updated"

This reverts commit 2e8ac3eefca27690c39fc29c17b7a5d31ad9c4ff.
---
 Project.toml     | 2 +-
 docs/make.jl     | 1 -
 test/runtests.jl | 1 -
 3 files changed, 1 insertion(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9ab64e9..c9f3c07 100644
--- a/Project.toml
+++ b/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.2.0"
+version = "0.1.4"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
diff --git a/docs/make.jl b/docs/make.jl
index 81305e5..10a5400 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,4 +1,3 @@
-using Revise
 using EllipseSampling
 using Documenter, DocumenterCitations
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a0f6ab9..2bdb38b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,3 @@
-using Revise
 using EllipseSampling
 using Test
 using LinearAlgebra

From 25ffce028158e9b1708220e73dca36f9728fedf6 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Tue, 23 Jan 2024 10:33:55 +1300
Subject: [PATCH 6/7] Update to v0.2.0 without adding `using Revise` lines in
 `runtests.jl` and `make.jl`

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index c9f3c07..9ab64e9 100644
--- a/Project.toml
+++ b/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.4"
+version = "0.2.0"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"

From 23120022f34e8f11fa9f1e30733b8ef0dc88e1b9 Mon Sep 17 00:00:00 2001
From: JoelTrent <79883375+JoelTrent@users.noreply.github.com>
Date: Tue, 23 Jan 2024 10:38:51 +1300
Subject: [PATCH 7/7] Add tests for ellipse parameter calculation at a higher
 degree of freedom

---
 test/runtests.jl | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 2bdb38b..8f5d505 100644
--- a/test/runtests.jl
+++ b/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