virgil.priors
Priors that are invariant measures but not plain numpyro distributions: the
isotropic orientation priors, uniform in cos i (orbit inclination) and in
sin(lat) (a position on a sphere). Locations take Uniform, scales take
LogUniform, and longitudes and nodes are uniform in their angles. See
Which priors? for the full table and the reasons.
hierarchical_scales draws a set of scales, such as one error scale per
epoch, from a log-normal population with log-uniform hyperpriors, for use as
noise= terms tied to the parameters.
Priors that are invariant measures but not plain numpyro priors.
virgil's default priors are the Jeffreys priors under the groups that act on
each parameter (see the design notes): uniform for locations, log-uniform
(numpyro.distributions.LogUniform) for scales, and the Haar measure for
orientations. The Haar measure on rotations is uniform in the direction of
an axis on the sphere, which is not uniform in the angles virgil uses:
- an orbit's inclination
inc(degrees, 0 to 180) has density ∝ sin i, so that cos i is uniform; useIsotropicInclination; - a position on a sphere, such as a spot's latitude (radians), has density
∝ cos(lat), so that sin(lat) is uniform; use
IsotropicLatitude.
The longitude of the node, argument of periastron and spot longitude are
uniform in their angles, so numpyro.distributions.Uniform is right for
them.
hierarchical_scales gives a set of
scale factors (say one error scale per epoch) drawn from a log-normal
population whose median and spread are themselves sampled, with
log-uniform hyperpriors.
Both classes are numpyro distributions, so they work as entries of the
priors of fit and
numpyro_model: their support is an
interval, which both map to unconstrained coordinates with numpyro's
biject_to. fit optimizes them in their flat coordinate (cos i or
sin(lat), through flat_coordinate()), where the prior is constant, so
Levenberg–Marquardt applies and is fit's automatic choice.
IsotropicInclination
Bases: _InverseCDFPrior
The inclination of an isotropically oriented axis: density ∝ sin i.
The Haar measure on rotations gives an orbit normal (or a spin axis)
that is uniform on the sphere, so its inclination i has density
sin(i) / (cos(low) - cos(high)) on [low, high], that is, cos i
is uniform. This is the default prior on inc; Uniform(0, 180)
is not isotropic (it favours face-on orbits).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
low
|
float
|
Bounds in degrees, the unit of |
0.0
|
high
|
float
|
Bounds in degrees, the unit of |
0.0
|
Notes
Use IsotropicInclination() (0 to 180) whenever the data distinguish
i from 180 - i: astrometry or interferometric closure phases
give the sense of motion on the sky, so the orbit is prograde
(i < 90) or retrograde (i > 90). Use
IsotropicInclination(0, 90) when only |cos i| is identifiable,
for example with radial velocities alone, or with the sky-projected
image of a rotating star: the data cannot tell i from 180 - i,
and the density is the same shape on half the range.
Examples:
>>> import jax
>>> from virgil.priors import IsotropicInclination
>>> prior = IsotropicInclination(0.0, 90.0)
>>> samples = prior.sample(jax.random.key(0), (1000,))
>>> bool(((samples >= 0.0) & (samples <= 90.0)).all())
True
mean
property
__init__(low=0.0, high=180.0, *, validate_args=None)
icdf(q)
Quantile function: cos i is uniform, inverted near the nearer pole.
cdf(value)
log_prob(value)
IsotropicLatitude
Bases: _InverseCDFPrior
The latitude of an isotropically placed point: density ∝ cos(lat).
A point uniform on the sphere has sin(lat) uniform, so its latitude has
density cos(lat) / (sin(high) - sin(low)). Use it for the position
of a star spot (harmonix's lat argument), together with
Uniform(-pi, pi) for the longitude.
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
low
|
float
|
Bounds in radians, the unit harmonix's |
-pi / 2
|
high
|
float
|
Bounds in radians, the unit harmonix's |
-pi / 2
|
Examples:
>>> import jax
>>> from virgil.priors import IsotropicLatitude
>>> samples = IsotropicLatitude().sample(jax.random.key(0), (1000,))
>>> bool((abs(samples) <= 3.1416 / 2).all())
True
mean
property
__init__(low=-np.pi / 2, high=np.pi / 2, *, validate_args=None)
icdf(q)
cdf(value)
log_prob(value)
PopulationScale
dataclass
Member index of a log-normal population of scales (see
hierarchical_scales).
Called with the dict of sampled values, it returns the member's scale
s_k. With centred=True, log s_k is sampled itself (at
"<name>_log", under a flat prior) and :meth:log_prior is its
population density, log N(log s_k | log median, spread), which
numpyro_model and
fit add once for each member used as a
noise term. Otherwise s_k = median exp(spread z_k) with z_k ~ N(0,
1) sampled at "<name>_z", and the log prior is zero.
It is a frozen dataclass, so equal members compare (and hash) equal: repeated fits do not recompile, and a member used for two terms adds its density once.
name
instance-attribute
index
instance-attribute
centred = True
class-attribute
instance-attribute
__init__(name, index, centred=True)
__call__(values)
log_prior(values)
The member's population log density (zero when non-centred, where the standard normal prior on z_k carries it).
hierarchical_scales(name, n, median=None, spread=None, centred=True)
n positive scales drawn from a log-normal population.
log s_k is normal about log median with standard deviation
spread, and median and spread are sampled too. They are
scales, so their default (Jeffreys) priors are log-uniform.
Use it for calibration nuisances that differ from epoch to epoch but
come from one instrument, such as one closure-phase error scale per
night: the population pulls poorly constrained epochs towards the
typical value, and its median and spread say how well the stated
errors describe the instrument. Pass the scales as noise terms
tied to the parameters, for
numpyro_model or
fit. Use every member: an unused member of
the centred form has a flat, improper prior.
Centred or not. The centred form (default) samples u_k = log s_k and adds the population density N(u_k | log median, spread). It suits members that the data measure well, such as error scales of epochs with tens of closure phases or more: u_k is then nearly independent of the hyperparameters. The non-centred form samples z_k ~ N(0, 1) with s_k = median exp(spread z_k), and suits members the data barely constrain; with well-measured members it makes a curved ridge, z_k ∝ 1/spread, on which NUTS diverges (in a test with eight members measured to 7%, 10 divergences in 4000 draws against none centred).
Parameters:
| Name | Type | Description | Default |
|---|---|---|---|
name
|
str
|
Prefix of the new parameters, |
required |
n
|
int
|
Number of scales, e.g. the number of epochs. |
required |
median
|
Distribution
|
Prior on the population median (default |
None
|
spread
|
Distribution
|
Prior on the standard deviation of log s (default
|
None
|
centred
|
bool
|
The parameterization (see above). |
True
|
Returns:
| Name | Type | Description |
|---|---|---|
priors |
dict
|
The three priors, to merge into a |
scales |
list of PopulationScale
|
One callable per member, each mapping the sampled values to s_k. |
Examples:
>>> import numpyro.distributions as dist
>>> from virgil.priors import hierarchical_scales
>>> priors, scales = hierarchical_scales("cp_scale", 3)
>>> sorted(priors)
['cp_scale_log', 'cp_scale_median', 'cp_scale_spread']
>>> noise = [{"phi_scale": s} for s in scales]
>>> values = {"cp_scale_median": 2.0, "cp_scale_spread": 0.5,
... "cp_scale_log": np.log(np.array([1.0, 2.0, 4.0]))}
>>> round(float(scales[1](values)), 6)
2.0