virgil.angles
Angles sampled as 2-D vectors, so that no prior has a wrap boundary. Put an
AngleVector in a priors dict under an angle's path (in degrees, such as a
binary's "pa" or an orbit's node), and fit, gauss_newton_mass and
numpyro_model sample a vector v = r(cos θ, sin θ) at the site
"<path>_vec" instead of the angle. The model receives θ in degrees, in
[0°, 360°): numpyro_model records it as the deterministic site "<path>",
and fit reports both.
from virgil.angles import AngleVector
priors = {"pa": AngleVector()} # uniform, no wall at 0°/360°
priors = {"pa": AngleVector(350.0, 20.0)} # von Mises, mean 350°, κ = 20
priors = {"node": AngleVector(100.0, 4.0, axial=True)} # known modulo 180°
The radius has a ring prior, ∝ exp(−(r − 1)²/2s²), whose mode is the unit
circle, so that maximum a posteriori fits stay away from the origin, where
θ is undefined. A von Mises prior enters as the chord √κ (v̂ − m̂), the same
form as virgil's unprojected phase residuals (2 sin(Δ/2)/σ, with κ = 1/σ²),
so every term has a least-squares form and Levenberg–Marquardt takes it. The
density is normalized in the plane (with the von Mises normalizer on the
circle, through i0e), so evidences stay normalized. The default, a uniform
angle, is the invariant prior; a von Mises prior is strong information from
an external measurement.
For orbits, orientation_priors samples 2Ω and ϖ = Ω + ω for
position-only fits (so the node ambiguity is one point, not two modes), or Ω
and ϖ when RVs fix the node.
Credit. The construction follows Octofitter's UniformCircular
(Thompson et al. 2023, AJ 166, 164) and exoplanet's Angle
(Foreman-Mackey et al. 2021, JOSS 6, 3285), which sample v ~ N(0, I). The
ring and the von Mises chords are virgil's. No code is taken from either.
Angles sampled as 2-D vectors, so that no prior has a wrap boundary.
A prior on an angle θ (degrees) with support on an interval, such as
Uniform(0, 360), numpyro's VonMises or
AxialVonMises, reaches samplers and
fit through a bijection onto the real line, so an
angle whose posterior straddles the end of the interval meets a wall there.
AngleVector removes the wall: put it in a
priors dict under the angle's path, and a vector v = r(cos θ, sin θ) in
ℝ² is sampled at the site "<path>_vec" instead. The angle itself, in
degrees in [0°, 360°), is what the model receives, and what
numpyro_model records as a deterministic
site "<path>" and fit reports under "<path>" (with the vector
under "<path>_vec").
The radius r is a nuisance with a ring prior, ∝ exp(-(r - 1)²/2s²), so the
mode is the unit circle and a maximum a posteriori fit keeps away from the
origin, where θ is undefined. The prior is rotationally symmetric, so θ is
exactly uniform unless a von Mises (or axial von Mises) prior is given;
that one enters as a chord, √κ (v̂ - m̂), the same form as virgil's
unprojected phase residuals (see
whitened_residuals). Every term has
a least-squares form, so Levenberg–Marquardt and
gauss_newton_mass take it unchanged,
and the density is normalized in ℝ², so evidences stay normalized.
The construction follows Octofitter's UniformCircular (Thompson et al.
2023, AJ 166, 164) and exoplanet's Angle (Foreman-Mackey et al. 2021,
JOSS 6, 3285), which sample v ~ N(0, I); the ring and the chord priors are
virgil's.
AngleVector
Bases: Distribution
A prior on an angle (degrees), sampled as a 2-D vector.
Put it in a priors dict, keyed by the angle's path, for
fit,
gauss_newton_mass or
numpyro_model. They sample the
vector v at the site "<path>_vec" and give the model its direction
θ = atan2(v₂, v₁) in degrees, in [0°, 360°). A fit starts from the
unit vector of the angle's starting value (in degrees, from init or
the template), or from init["<path>_vec"].
The density of v = r(cos θ, sin θ) is
p(v) = exp(-(r - 1)²/2s²) / Z_s × p(θ),
with Z_s = ∫₀^∞ r exp(-(r - 1)²/2s²) dr, so that r and θ are independent and θ has exactly the density p(θ) (per radian):
- uniform, 1/2π, by default;
- a von Mises, exp(κ(cos(θ - μ) - 1)) / (2π i0e(κ)), given
meanandkappa; - with
axial=True, the von Mises of 2θ, so that θ and θ + 180° are equally likely (asAxialVonMises).
Its least-squares residuals (residuals) are the ring, (r - 1)/s,
and for a von Mises the chord √κ (v̂ - m̂), with v̂ = v/r and m̂ =
(cos μ, sin μ): half its square is κ(1 - cos(θ - μ)), the von Mises
exponent, just as virgil's phase residuals 2 sin(Δ/2)/σ are chords with
κ = 1/σ². The axial chord is the same on the doubled direction,
(cos 2θ, sin 2θ) = (x² - y², 2xy)/r².
The default, uniform θ, is the invariant prior for an angle. A von Mises prior is strong information, from an external measurement.
The construction follows Octofitter's UniformCircular (Thompson et
al. 2023, AJ 166, 164) and exoplanet's Angle (Foreman-Mackey et al.
2021, JOSS 6, 3285), whose v ~ N(0, I) peaks at the origin, where θ is
undefined: a maximum a posteriori fit would drive r → 0. The ring's mode
is the unit circle instead, and every term has a least-squares form.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
mean
|
float
|
Mean angle μ (degrees) of a von Mises prior; with |
None
|
kappa
|
float
|
Concentration κ of the von Mises prior (of 2θ when |
None
|
axial
|
bool
|
Whether the von Mises prior is on 2θ (an axis known modulo 180°). |
False
|
ring_width
|
float
|
The ring's radial width s (default 0.25). It does not change the prior on θ. |
0.25
|
arg_constraints = {'ring_width': constraints.positive}
class-attribute
instance-attribute
pytree_data_fields = ('mean_deg', 'kappa', 'ring_width')
class-attribute
instance-attribute
pytree_aux_fields = ('axial',)
class-attribute
instance-attribute
support = constraints.real_vector
class-attribute
instance-attribute
mean_deg = None if mean is None else np.asarray(mean, float)
instance-attribute
kappa = None if kappa is None else np.asarray(kappa, float)
instance-attribute
ring_width = np.asarray(ring_width, float)
instance-attribute
axial = bool(axial)
instance-attribute
uniform
property
Whether the prior on the angle is uniform.
__init__(mean=None, kappa=None, *, axial=False, ring_width=0.25, validate_args=None)
residuals(vector)
Residuals whose half sum of squares is -log p(vector) + const.
The ring, (r - 1)/s, then, for a von Mises prior, the chord √κ (v̂ - m̂) (two more rows).
log_prob(value)
sample(key, sample_shape=())
vector_angle(vector)
The direction of vector (last axis of length 2), in degrees in
[0°, 360°), measured from its first component towards its second.
vector_site(path)
The numpyro site of the vector that carries the angle at path.
is_angle_vector(prior)
Whether prior is an AngleVector.