Single Stellar Populations

SSPModel(isofunc, data::AbstractMatrix{<:Number}, edges, mag_err_funcs, y_index, color_indices, imf,
         completeness_funcs, bias_funcs;
         logAge, MH, dmod, Av=0, binary_model::AbstractBinaryModel=NoBinaries(),
         binary_fraction=binary_system_fraction(binary_model), background_fraction=Uniform(0, 1),
         background=nothing, background_floor::Real=0.05, mask=falses(size(data)),
         mean_mass::Number=mean(imf))

Model of the observed Hess diagram data as a single simple stellar population (SSP) of one age and metallicity plus a background, such as a star cluster and its surrounding field. Fit with fit_ssp and sample with sample_ssp.

Arguments

• isofunc(logAge, MH, Av) returns the isochrone (m_ini, mags) for logarithmic age logAge, metallicity MH, and V-band extinction Av, with m_ini and mags as in partial_cmd_smooth. It is called from multiple threads, so it must be thread-safe, and it must support every value allowed by the priors on these parameters.
• data is the observed Hess diagram with bin edges edges, e.g., from bin_cmd.
• mag_err_funcs, y_index, color_indices, imf, completeness_funcs, and bias_funcs are as in partial_cmd_smooth, for the magnitudes returned by isofunc.

Parameters

The keyword arguments logAge, MH, dmod (distance modulus), Av, binary_fraction, and background_fraction are the parameters of the model. Each is either a Real, which fixes the parameter to that value, or a Distributions.UnivariateDistribution, which makes it a free parameter with that distribution as its prior. Free parameters are ordered as in model.free. binary_fraction is the binary system fraction of binary_model (see with_fraction), and background_fraction is the fraction of the observed stars in the unmasked bins of data that belong to the background.

Keyword Arguments

• binary_model::AbstractBinaryModel is the binary model, whose binary fraction is replaced by binary_fraction.
• background is the shape of the background, a matrix the same size as data (e.g., the Hess diagram of a nearby field), or nothing for a uniform background.
• background_floor is the fraction of the background spread uniformly over the unmasked bins, so that stars in bins where both background and the SSP are zero can be explained by the background rather than distorting the SSP parameters.
• mask is true for bins excluded from the fit (see hess_mask).
• mean_mass is as in partial_cmd_smooth.

Likelihood

The expected number of stars in bin $i$ is $A \, p_i$ with

\[p_i = (1 - g) \frac{c_i}{\sum_j c_j} + g \, u_i\]

where $c_i$ is the SSP Hess diagram per solar mass formed (see ssp_components), $u_i$ is the background shape normalized to sum to 1, $g$ is background_fraction, and sums run over the unmasked bins. With a scale-invariant prior $p(A) \propto 1/A$ on the total expected number of stars $A$, marginalizing over $A$ gives the multinomial log-likelihood $\sum_i n_i \log p_i$ for the observed counts $n_i$, so the fit depends only on the shape of the Hess diagram. The posterior of $A$ is then $\text{Gamma}(N, 1)$ for $N = \sum_i n_i$ observed stars, from which sample_ssp derives the stellar mass formed in the SSP. This model implements the LogDensityProblems interface for use with other samplers.

fit_ssp(model::SSPModel; x0=nothing, ncandidates::Integer=1000, nstarts::Integer=16,
        options::Optim.Options=Optim.Options(iterations=1000, g_tol=0.01), restarts::Integer=3,
        restart_tol::Real=1, rng::AbstractRNG=default_rng())::SSPFit

Finds the maximum a posteriori (MAP) parameters of model by local optimization with the Nelder-Mead method from one or more starting points. The optimization is done in unconstrained coordinates, so trial points stay within the support of the priors.

If x0, a vector of the free parameters in the order of model.free, is given, it is the only starting point. Otherwise, the log posterior density is evaluated at ncandidates candidate starting points spread evenly over the free logAge, MH, and dmod, each between the 0.02 and 0.98 quantiles of its prior (0.001 and 0.999 for a side on which the prior is unbounded), by the quasirandom R2 sequence [10] with a random offset drawn from rng. The other free parameters are at their prior medians, except that free binary and background fractions are maximized at each candidate (which does not require recomputing the templates). The local optimization starts from every candidate that is the best in its neighborhood among the best 20% of the candidates, by multi-level single linkage [11], which finds one starting point for each separate peak the candidates resolve. These run in parallel over threads. The best of their solutions is the MAP, and the solution of every local optimization is kept in starts (one per starting point, so two starting points can reach the same peak). Multiple starting points guard against converging to a secondary maximum, which the degeneracies between age, metallicity, and distance can produce. If there are more than nstarts starting points, only the best nstarts are used, with a warning. The extinction is not searched but starts at its prior median, so give x0 if its prior is much wider than about 0.2 mag. The result depends on rng but not on the number of threads.

Because Nelder-Mead can stall before reaching the maximum, each local optimization is restarted from its solution with a fresh simplex, up to restarts times, until the log posterior density improves by less than restart_tol. options are passed to Optim.optimize; by default, each run stops when the standard deviation of the negative log density over the simplex falls below g_tol = 0.01. The default tolerances are loose because the MAP only needs to be within the region of the posterior that sample_ssp explores; posterior samples typically lie a few units of log density below the maximum. Returns an SSPFit; sample the posterior with sample_ssp, which uses every local solution in starts to find the modes of the posterior.

SSPFit

Result of fit_ssp, with fields

• model::SSPModel, the fitted model,
• map::NamedTuple, the value of every parameter at the maximum a posteriori (MAP) solution,
• θ::Vector{Float64}, the free parameters at the MAP, in the order of model.free,
• logdensity::Float64, the log posterior density at the MAP (up to a constant),
• candidates, nothing or a NamedTuple with the free parameters at each candidate starting point of fit_ssp, with the free binary and background fractions maximized (θ), and the log posterior density there (logdensity),
• starts, a vector with one NamedTuple per local optimization, with its starting point θ0, solution θ, and log posterior density logdensity,
• optim, the Optim.jl results of the Nelder-Mead runs of the best local optimization.

sample_ssp(fit::SSPFit; pilot_steps::Integer=300, pilot_walkers::Integer=max(16, 4 * length(fit.θ)),
           is_rounds::Integer=3, is_samples::Integer=2000, nsamples::Integer=is_rounds * is_samples,
           mode_cut::Real=20, radius::Real=0.5, rng::AbstractRNG=default_rng())::SSPSamples

Samples the posterior of the model fit by fit_ssp, including every mode found by its local optimizations, and returns an SSPSamples with equal-weight samples, the posterior probability and evidence of each mode, and convergence diagnostics.

1. Modes. The local solutions in fit.starts within mode_cut of the best log posterior density are grouped into modes. From each, an ensemble of pilot_walkers walkers (an even number, at least the number of free parameters plus 2) of the affine-invariant sampler of KissMCMC.jl, started in a ball of radius radius times the posterior width along each parameter, takes pilot_steps steps; the mean and covariance of the second half describe the mode. Two solutions are the same mode if either lies within the 99.9% chi-squared distance (for the number of free parameters) of the other's pilot, so solutions that belong to an existing mode usually need no pilot of their own; when one does, its pilot is pooled with the mode's, so that solutions spread along an elongated mode together describe it.
2. Adaptive importance sampling. The posterior is sampled with is_rounds rounds of is_samples draws each from a mixture of multivariate $t$ distributions, one per mode, starting from the pilots' means and covariances. After each round, each component is refitted to the importance-weighted draws it is responsible for. At least 3 rounds are required, since fewer often fail to adapt the proposal from the pilots. The draws of all rounds are combined with the weights of adaptive multiple importance sampling (AMIS), which give each mode's posterior probability and evidence directly.
3. Samples. nsamples equal-weight samples are drawn from the weighted draws by systematic resampling, so they contain repeated draws. The total expected number of stars is drawn from its posterior for each sample (see SSPModel) to give the stellar mass formed and the expected number of background stars.

The importance weights also tell whether the settings were adequate. A warning is issued, and recorded in the result, if the last round's proposal was still adapting (a proposal shift above 0.5 nats or a Pareto $\hat k$ of the weights above 0.7), in which case is_rounds or pilot_steps should be increased, or if the weights have fewer than 1000 effective samples, in which case is_samples should be increased. The evidence and mode probabilities are relative to the modes that fit_ssp found; a mode that none of its local optimizations reached is not represented.

Evaluations of the log posterior density are parallelized over threads. Results are reproducible for a given rng and number of threads.

SSPSamples

Result of sample_ssp, with fields

• fit::SSPFit, the fit that was sampled,
• samples::NamedTuple, equal-weight posterior samples as one vector per column: the free parameters, named as in fit.model.free; mass, the stellar mass formed in the SSP in solar masses; background_stars, the expected number of background stars in the unmasked bins; lp, the log posterior density (up to a constant); and mode, the index into modes of the mode each sample belongs to,
• logZ::Float64, the logarithm of the evidence, the integral of the posterior density over the parameters, up to the same constant as lp,
• ess::Float64, the effective sample size of the importance weights; the samples contain repeated draws, and ess is the number of independent samples they are worth,
• modes, a vector with one NamedTuple per mode: the local solution of fit_ssp it was found from (θ and logdensity), its posterior probability (weight), its log evidence (logZ), and the effective sample size of its importance weights (ess),
• diagnostics, a NamedTuple with the Pareto $\hat k$ of the importance weights of the last round (khat), the proposal shift of each mode in the last round (shift, in nats), and the log evidence and fraction of effective draws of each round (rounds),
• warnings::Vector{String}, the warnings issued by sample_ssp.

ssp_hess(model::SSPModel, p::NamedTuple)
ssp_hess(fit::SSPFit)

Returns the expected Hess diagram of model for the parameters p (a value for every parameter, e.g., fit.map from fit_ssp) as a matrix the same size as the observed Hess diagram. It is normalized to the number of observed stars in the unmasked bins, which is the posterior mean of the total expected number of stars. Masked bins are included, so the model can be compared to the data inside them.

ssp_mass(model::SSPModel, p::NamedTuple)
ssp_mass(fit::SSPFit)

Returns the stellar mass formed in the SSP of model, in solar masses, for the parameters p (a value for every parameter, e.g., fit.map from fit_ssp), with the total expected number of stars equal to the number of observed stars in the unmasked bins, as in ssp_hess. sample_ssp instead draws the total expected number of stars from its posterior for each sample, so its mass samples also include the Poisson uncertainty in the number of stars.

(single, binary) = ssp_components(m_ini::AbstractVector{<:Number},
                                  mags::AbstractVector{<:AbstractVector{<:Number}},
                                  mag_err_funcs, y_index, color_indices, imf,
                                  completeness_funcs=[one for i in mags],
                                  bias_funcs=[zero for i in mags];
                                  binary_model::AbstractBinaryModel=NoBinaries(),
                                  dmod::Number=0, normalize_value::Number=1,
                                  mean_mass::Number=mean(imf),
                                  edges=nothing, xlim=nothing, ylim=nothing, nbins=nothing,
                                  xwidth=nothing, ywidth=nothing)

Returns the Hess diagram of the simple stellar population as two components, both StatsBase.Histograms normalized to normalize_value solar masses formed: single, for a population of only single stars, and binary, for a population of only binary systems under binary_model (i.e., with_fraction(binary_model, 1)). binary is nothing for binary_model = NoBinaries(). Arguments and keyword arguments are as in partial_cmd_smooth.

The Hess diagram for a binary system fraction b is (1 - φ) .* single.weights .+ φ .* binary.weights with φ = binary_mass_fraction(with_fraction(binary_model, b), imf), which is equal to the result of partial_cmd_smooth for with_fraction(binary_model, b). These components therefore give the Hess diagram for any binary fraction without recomputing them.

(single, binary) = ssp_components(m_ini::AbstractVector{<:Number},
                                  mags::AbstractVector{<:AbstractVector{<:AbstractVector{<:Number}}},
                                  knots::AbstractVector{<:Number},
                                  extinction::Distributions.ContinuousUnivariateDistribution,
                                  mag_err_funcs, y_index, color_indices, imf,
                                  completeness_funcs=[one for i in first(mags)],
                                  bias_funcs=[zero for i in first(mags)]; kws...)

Method of ssp_components for a distribution of extinctions; arguments are as for the corresponding method of partial_cmd_smooth.

adaptive_is(logd, μs, Σs; rounds, n, ν=5, inflate=1.5, minweight=0.1, rng)

Adaptive importance sampling of the log density logd with the mixture population Monte Carlo (M-PMC) method of Cappé et al. [12], whose proposal is a mixture of multivariate $t$ distributions that is refitted to the importance-weighted draws after each of rounds rounds. logd(x) returns the log density (normalized or not) at the parameter vector x and a Float64 blob that is stored with each draw.

The mixture has one component per mode, with ν > 2 degrees of freedom; component $k$ starts with weight $\alpha_k = 1/K$, where $K$ is length(μs), mean μs[k], and covariance inflate^2 * Σs[k]. Each round

1. draws n points $x_i$ from the mixture with rng;
2. evaluates logd at the draws in parallel over threads;
3. weights each draw by $w_i \propto \pi(x_i) / q(x_i)$, where $\pi$ is the target density exp(logd) and $q(x) = \sum_k \alpha_k \, t_k(x)$ is the density of the whole mixture;
4. refits the mixture as in the M-PMC update of Cappé et al. [12]: the draws are split among the components by their responsibilities $\rho_{ik} = \alpha_k t_k(x_i) / q(x_i)$, each $\alpha_k$ becomes $\sum_i w_i \rho_{ik}$ (with normalized $w_i$), and each mean and covariance becomes the mean and covariance of the draws weighted by $w_i \rho_{ik}$.

The update differs from that of Cappé et al. [12] in ways that keep the proposal robust when it starts far from the target. The means and covariances are refitted as weighted moments, rather than with their update for $t$ mixtures, and the covariances are again multiplied by inflate^2, so the proposal keeps heavier tails than the target, which bounds the importance weights. Each $\alpha_k$ is kept at least minweight (before renormalizing), so that no mode is starved of draws. A component whose weighted draws are worth fewer than twice the number of parameters in effective samples keeps its previous mean and covariance.

Returns a vector with one NamedTuple per round with fields Y (a matrix with the draws as columns), ld and blob (the values returned by logd at the draws), lw (the log importance weights, $\log \pi(x_i) - \log q(x_i)$), and the proposal of that round: α (the component weights), μs (the means), and Λs (the scale matrices, which are the covariances times $(\nu - 2) / \nu$). amis_combine combines the draws of all rounds.

amis_combine(rounds, ν)

Combines the draws of all rounds returned by adaptive_is, whose mixture components have ν degrees of freedom, with the weights of adaptive multiple importance sampling (AMIS) [13]. Rather than by its own round's proposal, each draw $x_i$ is weighted by the mixture of the proposals $q_r$ of all rounds $r$, each weighted by its share of the $N$ draws, $N_r / N$,

\[w_i = \frac{\pi(x_i)}{\sum_r (N_r / N) \, q_r(x_i)}.\]

The early rounds, whose proposals were still adapting, then add to the sample instead of being discarded, and a round whose proposal missed part of the target is covered by the others. As Cornuet et al. [13] discuss, because each proposal depends on the earlier draws, the convergence of AMIS estimators does not follow from that of standard importance sampling.

Returns a NamedTuple with fields X (a matrix with the draws of all rounds as columns), ld and blob (the values returned by logd at the draws), lw (the log AMIS weights), and mode (for each draw, the component that contributes most to the combined mixture density at the draw, which assigns the draw to a mode).

pareto_khat(lw)

Returns the Pareto $\hat k$ diagnostic of Vehtari et al. [14] for the importance weights with logarithms lw: the shape of a generalized Pareto distribution fitted by gpd_khat to the largest $M = \min(S/5, 3\sqrt{S})$ of the $S$ finite weights, as exceedances over the next largest weight. It measures how heavy the tail of the weights is, which sets how reliable estimates from them are: they are reliable for $\hat k$ below 0.7, the threshold used by sample_ssp, and above it a few draws carry most of the weight and the estimates can be badly biased even when the effective sample size looks adequate. Returns Inf if there are too few finite weights to fit (fewer than 5 in the tail) or if the largest weight is so much larger than the rest of the tail that they underflow to 0.

gpd_khat(x)

Estimates the shape $k$ of a generalized Pareto distribution from the nonnegative exceedances x over a threshold, with the method of Zhang and Stephens [15], which reparameterizes the distribution by the ratio $\theta$ of its shape to its scale (up to sign), estimates $\theta$ as the mean over a grid of $m = 30 + \lfloor \sqrt{n} \rfloor$ values, for $n$ exceedances, weighted by their profile likelihoods (an approximate posterior mean), and takes $k$ as its maximum-likelihood value given $\theta$. As in the implementation of Pareto smoothed importance sampling (PSIS) by Vehtari et al. [14], $k$ is then shrunk toward 0.5 by a weakly informative prior worth 10 observations, which stabilizes the estimate for small $n$. Larger $k$ means a heavier tail; the distribution has a finite mean for $k < 1$ and a finite variance for $k < 1/2$.