Merge pull request #34 from JoelTrent/dev

Update to version 0.1.3

Project.toml

@@ -1,11 +1,12 @@
 name = "EllipseSampling"
 uuid = "7d220c50-6298-48cd-9710-1d1a8ef13806"
 authors = ["JoelTrent <79883375+JoelTrent@users.noreply.github.com> and contributors"]
-version = "0.1.2"
+version = "0.1.3"
 
 [deps]
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Elliptic = "b305315f-e792-5b7a-8f41-49f472929428"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SnoopPrecompile = "66db9d55-30c0-4569-8b51-7e840670fc0c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

docs/Manifest.toml

@@ -26,7 +26,7 @@ version = "0.9.3"
 deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
 git-tree-sha1 = "58fea7c536acd71f3eef6be3b21c0df5f3df88fd"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.24"
+version = "1.0.0"
 
 [[deps.EllipseSampling]]
 path = ".."

docs/Project.toml

@@ -1,3 +1,4 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 EllipseSampling = "7d220c50-6298-48cd-9710-1d1a8ef13806"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"

docs/make.jl

@@ -1,12 +1,16 @@
 using EllipseSampling
-using Documenter
+using Documenter, DocumenterCitations
 
 DocMeta.setdocmeta!(EllipseSampling, :DocTestSetup, :(using EllipseSampling); recursive=true)
 
+bib = CitationBibliography(
+    joinpath(@__DIR__, "src", "refs.bib");
+    style=:numeric
+)
+
 makedocs(;
     modules=[EllipseSampling],
     authors="JoelTrent <79883375+JoelTrent@users.noreply.github.com> and contributors",
-    repo="https://github.com/JoelTrent/EllipseSampling.jl/blob/{commit}{path}#{line}",
     sitename="EllipseSampling.jl",
     format=Documenter.HTML(;
         prettyurls=get(ENV, "CI", "false") == "true",
@@ -18,8 +22,10 @@ makedocs(;
         "Home" => "index.md",
         "Quick Start" => "quick_start.md",
         "User Interface" => "user_interface.md",
-        "Internal Library" => "internal_library.md"
+        "Internal Library" => "internal_library.md",
+        "References" => "references.md"
     ],
+    plugins=[bib],
 )
 
 deploydocs(;

docs/src/references.md

@@ -0,0 +1,4 @@
+# References
+
+```@bibliography
+```

docs/src/refs.bib

@@ -0,0 +1,24 @@
+@book{pawitanall2001,
+	title = {In All Likelihood: Statistical Modelling and Inference Using Likelihood},
+	isbn = {978-0-19-850765-9},
+	shorttitle = {In All Likelihood},
+	pagetotal = {528},
+	publisher = {Oxford University Press},
+	author = {Pawitan, Yudi},
+	year = {2001},
+}
+
+
+@article{friendlyelliptical2013,
+	title = {Elliptical Insights: Understanding Statistical Methods through Elliptical Geometry},
+	volume = {28},
+	issn = {0883-4237},
+	url = {https://www.jstor.org/stable/43288410},
+	shorttitle = {Elliptical Insights},
+	pages = {1--39},
+	number = {1},
+	journal = {Statistical Science},
+	author = {Friendly, Michael and Monette, Georges and Fox, John},
+	year = {2013},
+	note = {Publisher: Institute of Mathematical Statistics},
+}

src/EllipseSampling.jl

@@ -1,7 +1,7 @@
 module EllipseSampling
 
 import Roots, Elliptic, Distributions
-using StaticArrays
+using StaticArrays, LinearAlgebra
 
 export construct_ellipse,
     generate_N_equally_spaced_points, generate_perimeter_point,

src/mathematic_relationships.jl

@@ -11,7 +11,6 @@ Computes `m`, the eccentricity of the ellipse squared, ``m=e^2``, where ``m=1-\\
 This relationship between m and e is seen by considering the equation for e: ``e=\\frac{c}{a}``, where ``c^2=\\big(a^2-b^2\\big)`` and ``a>b``.
 
 Replacing c in the equation for e: 
-    
 ```math
 e=\\frac{\\sqrt{a^2-b^2}}{a}
 ```

src/matrix_to_parameters.jl

@@ -5,15 +5,42 @@
 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.
 
 # Arguments
-- `Γ`: a square matrix (2D) which is the inverse of the Hessian of a log-likelihood function at its maximum likelihood estimate.
+- `Γ`: 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.
 
 # 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 2π radians that the ellipse has been rotated by. `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`.
+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`.
 
-References and equation(s) to come.
+## Observed Fisher Information Matrix and Approximation of the Log-Likelihood Function
+
+We can approximate a log-likelihood function from multiple parameters, ``\\theta``, by considering the observed Fisher information matrix (FIM). The observed FIM is
+a quadratic approximation of the curvature of the log-likelihood function at the maximum likelihood estimate, ``\\hat{\\theta}``. The observed FIM is the matrix of second derivatives (the Hessian) of the log-likelihood function evaluated at the MLE with elements [pawitanall2001](@cite): 
+```math
+H_{jk}(\\hat{\\theta}) \\equiv -\\frac{\\partial^2}{\\partial \\theta_j \\partial \\theta_k} \\ell (\\hat{\\theta} \\, ; \\, y_{1:I}^{\\textrm{o}} ).
+```
+
+This then allows us to define the following approximation of the normalised log-likelihood function
+using a second-order Taylor expansion at the MLE [pawitanall2001](@cite):
+```math
+    \\hat{\\ell} (\\theta \\, ; \\, y_{1:I}^{\\textrm{o}} ) \\approx -\\frac{1}{2} (\\theta-\\hat{\\theta})' H(\\hat{\\theta}) (\\theta-\\hat{\\theta}).
+```
+
+Similarly, for two parameters, ``\\psi``, from a larger number of parameters, we first invert the observed FIM, ``\\Gamma(\\hat{\\theta}) = H^{-1}(\\hat{\\theta})``, and then select just the rows and columns relating to the parameters of interest, before again inverting the matrix:
+```math
+    \\hat{\\ell}_p (\\psi \\, ; \\, y_{1:I}^{\\textrm{o}} )  \\approx -\\frac{1}{2} (\\psi-\\hat{\\psi})' ([e_j, e_k]' \\, \\Gamma(\\hat{\\theta}) \\, [e_j, e_k])^{-1} (\\psi-\\hat{\\psi}), \\hspace{0.2cm} \\theta_j \\cup \\theta_k = \\psi, 
+```
+where ``e_j`` and ``e_k`` are the ``j``th and ``k``th canonical vectors of ``\\mathbb{R}^{|\\theta|}``.
+
+## 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):
+```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)
@@ -24,37 +51,15 @@ function calculate_ellipse_parameters(Γ::Matrix{Float64}, ind1::Int, ind2::Int,
     (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)
-
-    α = atan(2*Hw[1,2]/(Hw[1,1]-Hw[2,2]))/2
-
-    # if α close to +/-0.25pi, +/-1.25pi, then switch to BigFloat precision
-    if isapprox(abs(rem(α/pi, 1)), 0.25, atol=1e-2)
-
-        # convert values to BigFloat for enhanced precision - required for correct results when α → 0.25pi or 1.25pi.
-        Hw = inv(BigFloat.(Γ[[ind1, ind2], [ind1, ind2]], RoundUp, precision=64)) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5) 
-
-        α = atan(2*Hw[1,2]/(Hw[1,1]-Hw[2,2]))/2
-        y_radius = sqrt( (cos(α)^4 - sin(α)^4) / (Hw[2,2]*(cos(α)^2) - Hw[1,1]*(sin(α)^2))  )
-        x_radius = sqrt( (cos(α)^2) / (Hw[1,1] - (sin(α)^2)/y_radius^2))
-
-        α, x_radius, y_radius = convert(Float64, α), convert(Float64, x_radius), convert(Float64, y_radius)
-    else
-        y_radius = sqrt( (cos(α)^4 - sin(α)^4) / (Hw[2,2]*(cos(α)^2) - Hw[1,1]*(sin(α)^2))  )
-        x_radius = sqrt( (cos(α)^2) / (Hw[1,1] - (sin(α)^2)/y_radius^2))
-    end
-
-    a = max(x_radius, y_radius)
-    b = min(x_radius, y_radius)
-
-    # a, b and α could also be calculated by finding the sqrt of the reciprocal of the eigenvalues of the normalised inv(Γ)
-    # using LinearAlgebra
-    # sqrt.(1.0 ./ eigvals(inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5)))
-
-    # this method for finding α assumes that a is the x axis
-    # eigs = eigvecs(inv(Γ[[2,3], [2,3]]) .* 0.5 ./ (quantile(Chisq(2), 0.95)*0.5))
-    # atan(eigs[2,1], eigs[1,1]) # if a is x axis
-    # atan(eigs[2,1], eigs[1,1]) + pi/2 # if a is y axis
-
-    return a, b, x_radius, y_radius, α 
+    Hw = inv(Γ[[ind1, ind2], [ind1, ind2]]) .* 0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level) * 0.5)
+    eigs = eigen(Hw)
+    a_eig, b_eig = 1.0 ./ sqrt.(eigs.values)
+
+    # α_eig constructed such that it is the angle in radians from the x axis to the major axis, 
+    # i.e. angle between the x axis and the largest eigenvector
+    α_eig = atan(eigs.vectors[2,1], eigs.vectors[1,1]) # value in interval [-pi, pi]
+    if α_eig < 0.0; α_eig += π end # value in interval [0, pi]
+
+    x_radius = a_eig; y_radius=b_eig
+    return a_eig, b_eig, x_radius, y_radius, α_eig
 end

test/runtests.jl

@@ -132,11 +132,10 @@ end
         @test isapprox_ellipsesampling(a_eig, a)
         @test isapprox_ellipsesampling(b_eig, b)
 
-        if x_radius > y_radius
-            @test isapprox_ellipsesampling(atan(eigvectors[2,1], eigvectors[1,1]), α)
-        else
-            @test isapprox_ellipsesampling(atan(eigvectors[2,1], eigvectors[1,1]) + 0.5*pi, α)
-        end
+        α_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 
@@ -151,10 +150,11 @@ end
         Hw = Hw_norm ./ (0.5 ./ (Distributions.quantile(Distributions.Chisq(2), confidence_level)*0.5))
         Γ = convert.(Float64, inv(BigFloat.(Hw, precision=64)))
 
-        a_calc, b_calc, _, _, _ = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level)
+        a_calc, b_calc, _, _, α_calc = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level)
 
         @test isapprox_ellipsesampling(a, a_calc)
         @test isapprox_ellipsesampling(b, b_calc)
+        @test isapprox_ellipsesampling(α, α_calc)
     end
 
 end