Minimum working example (for version 0.1.1):
using EllipseSampling
import Distributions
using Plots
gr()
N=8
a, b = 2.0, 1.0
α = 0.25*π
Cx, Cy= 2.0, 2.0
points=generate_N_equally_spaced_points(N, a, b, α, Cx, Cy, start_point_shift=0.0)
plt=scatter(points[1,:], points[2,:], label=nothing, mw=0, ms=8, aspect_ratio=:equal)
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(2), confidence_level)*0.5))
Γ = inv(Hw)
points2 = generate_N_equally_spaced_points(N, Γ, [Cx, Cy], 1, 2, confidence_level=confidence_level, start_point_shift=0.0)
scatter!(points2[1, :], points2[2, :], label=nothing, markeralpha=0.5, mw=0, ms=4)
display(plt)
If we modify the rotation slightly away from 0.25pi the problem goes away:
a, b = 2.0, 1.0
α = 0.251*π
Cx, Cy= 2.0, 2.0
points=generate_N_equally_spaced_points(N, a, b, α, Cx, Cy, start_point_shift=0.0)
plt=scatter(points[1,:], points[2,:], label=nothing, mw=0, ms=8, aspect_ratio=:equal)
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(2), confidence_level)*0.5))
Γ = inv(Hw)
points2 = generate_N_equally_spaced_points(N, Γ, [Cx, Cy], 1, 2, confidence_level=confidence_level, start_point_shift=0.0)
scatter!(points2[1, :], points2[2, :], label=nothing, markeralpha=0.5, mw=0, ms=4)
display(plt)
We can confirm this issue by inspecting the ellipse values returned by calculate_ellipse_parameters
:
# for α = 0.25*π
a_calc, b_calc, _, _, α_calc = EllipseSampling.calculate_ellipse_parameters(Γ, 1, 2, confidence_level)
println("a=", a_calc); println("b=", b_calc); println("α/pi=", α_calc/pi)
a=1.7994708216848752
b=1.0307764064044151
α=0.24999999999999992