Overview

AbstractMetallicityModel{T <: Real} is the abstract supertype for all hierarchical metallicity models. Abstract subtypes are AbstractAMR for age-metallicity relations and AbstractMZR for mass-metallicity relations.

fit_sfh(MH_model0::AbstractMetallicityModel,
        disp_model0::AbstractDispersionModel,
        models::AbstractMatrix{<:Number},
        data::AbstractVector{<:Number},
        logAge::AbstractVector{<:Number},
        metallicities::AbstractVector{<:Number};
        x0::AbstractVector{<:Number} = <...>,
        background::AbstractBackground = NoBackground(),
        iterations::Integer = 500,
        kws...)
fit_sfh(MH_model0::AbstractMetallicityModel,
        disp_model0::AbstractDispersionModel,
        models::AbstractVector{<:AbstractMatrix{<:Number}},
        data::AbstractMatrix{<:Number},
        logAge::AbstractVector{<:Number},
        metallicities::AbstractVector{<:Number};
        x0::AbstractVector{<:Number} = <...>,
        background::AbstractBackground = NoBackground(),
        iterations::Integer = 500,
        kws...)

Returns a CompositeBFGSResult instance that contains the maximum a posteriori (MAP) and maximum likelihood estimates (MLE) obtained from fitting the provided simple stellar population (SSP) templates models (with logarithmic ages logAge = log10(age [yr]) and metallicities metallicities) to the provided data. The metallicity evolution is modelled using the provided MH_model0, whose parameters can be free or fixed, with metallicity dispersion at fixed time modelled by disp_model0, whose parameters can be free or fixed.

This method is designed to work best with a grid of stellar models, defined by the outer product of N unique entries in logAge and M unique entries in metallicities. See the examples for more information on usage.

We provide several options for age-metallicity relations and mass-metallicity relations that can be used for MH_model0 and define APIs for users to create new models that will integrate with this function. Similar flexibility is allowed for the metallicity dispersion model disp_model0.

The primary method signature uses flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details, as well as stack_models that facilitates rearranging the models into this flattened format.

Arguments

• MH_model0 is an instance of AbstractMetallicityModel that defines how the average metallicity stars being formed in the population changes over time. The fittable parameters contained in this instance are used as the initial values to start the optimization.
• disp_model0 is an instance of AbstractDispersionModel that defines the distribution of metallicities of stars forming in a fixed time bin (i.e., the dispersion in metallicity around the mean at fixed time). The fittable parameters contained in this instance are used as the initial values to start the optimization.
• models are the template Hess diagrams for the SSPs that compose the observed Hess diagram.
• data is the Hess diagram for the observed data.
• logAge::AbstractVector{<:Number} is the vector containing the effective ages of the stellar populations used to create the templates in models, in units of log10(age [yr]). For example, if a population has an age of 1 Myr, its entry in logAge should be log10(10^6) = 6.0.
• metallicities::AbstractVector{<:Number} is the vector containing the effective metallicities of the stellar populations used to create the templates in models. These should be logarithmic abundances like [M/H] or [Fe/H]. There are some notes on the Wikipedia that might be useful.

Keyword Arguments

• x0 is the vector of initial guesses for the stellar mass coefficients per unique entry in logAge. We try to set reasonable defaults, but in most cases users should be calculating and passing this keyword argument. We provide StarFormationHistories.construct_x0_mdf to prepare x0 assuming a constant star formation rate and total stellar mass, which is typically a good initial guess.
• background is the background model, NoBackground (the default) or HessBackground, whose expected number of background stars is fit after the metallicity dispersion model parameters. If the HessBackground was constructed without an initial number of stars, the fit starts from 1% of the stars in data. With a background, x0 is renormalized so that the initial model, including the initial background, has as many stars as the data.
• iterations is the maximum number of BFGS iterations for each of the MAP and MLE optimizations; converged fits typically take 100-200.
• kws... are passed to Optim.Options and can be used to control tolerances for convergence.

Returns

• This function returns a CompositeBFGSResult that contains the output from both MLE and MAP optimizations, accessible via result.mle and result.map. These are each instances of BFGSResult. See the docs for these structs for more information.

CompositeBFGSResult(map::BFGSResult, mle::BFGSResult)

Type for containing the maximum a posteriori (MAP) AND maximum likelihood estimate (MLE) results from BFGS optimizations that use Optim.jl, which are individually accessible via the :mle and :map properties (i.e., for an instance of this type t, t.mle or getproperty(t, :mle) and t.map or getproperty(t, :map)).

Random samples can be drawn from an instance t as rand(t, N::Integer). This will return a size length(μ) x N matrix. This will use the MLE result for the best-fit values and the inverse Hessian approximation to the covariance matrix from the MAP result, which is more robust when best-fit values that are constrained to be positive approach 0.

Per-SSP coefficients can be calculated with calculate_coeffs(result::CompositeBFGSResult, logAge::AbstractVector{<:Number}, metallicities::AbstractVector{<:Number}), which uses the MLE result (see these docs).

BFGSResult(μ::AbstractVector{<:Number},
           σ::AbstractVector{<:Number},
           invH::AbstractMatrix{<:Number},
           result,
           MH_model::AbstractMetallicityModel,
           disp_model::AbstractDispersionModel,
           bg_model::AbstractBackground = NoBackground())

Type for containing the maximum likelihood estimate (MLE) or maximum a posteriori (MAP) results from BFGS optimizations that use Optim.jl. Fields are as follows:

• μ contains the final values of the fitting parameters. The mode and median methods will both return μ, but the mean of samples is not always equal to μ due to the variable transformations we perform.
• σ contains the standard errors estimated for the parameters and is returned by the std method.
• invH is the BFGS approximation to the inverse Hessian, which is an estimator for the covariance matrix of the parameters if the objective function is approximately Gaussian near the best-fit μ.
• result is the full result object returned by Optim.jl.
• MH_model is the best-fit metallicity model.
• disp_model is the best-fit metallicity dispersion model.
• bg_model is the best-fit background model (see HessBackground), whose parameters are the last elements of μ.

This type is implemented as a subtype of Distributions.Sampleable{Multivariate, Continuous} to enable sampling from an estimate of the likelihood / posterior distribution constructed from the invH. You can obtain N::Integer samples from the distribution with rand(R, N) where R is an instance of this type. This will return a size length(μ) x N matrix.

You can also directly obtain the per-SSP template coefficients ($r_{j,k}$ in the derivation) using the optimization results stored in a BFGSResult with calculate_coeffs.

See also

• CompositeBFGSResult is a type that contains two instances of BFGSResult, one for the MAP and one for the MLE.

sample_sfh(bfgs_result::CompositeBFGSResult, 
           models::AbstractMatrix{<:Number},
           data::AbstractVector{<:Number},
           logAge::AbstractVector{<:Number},
           metallicities::AbstractVector{<:Number},
           Nsteps::Integer;
           ϵ::Real = 0.05, # HMC step size
           reporter = DynamicHMC.ProgressMeterReport(),
           show_convergence::Bool = true,
           rng::AbstractRNG = default_rng())
sample_sfh(bfgs_result::CompositeBFGSResult, 
           models::AbstractVector{<:AbstractMatrix{<:Number}},
           data::AbstractMatrix{<:Number},
           logAge::AbstractVector{<:Number},
           metallicities::AbstractVector{<:Number},
           Nsteps::Integer;
           kws...)

Takes the SFH fitting result in bfgs_result and uses it to initialize the Hamiltonian Monte Carlo (HMC) sampler from DynamicHMC.jl to sample Nsteps independent draws from the posterior.

The primary method signature uses flattened formats for models and data. See the notes for the flattened call signature of StarFormationHistories.composite! for more details, as well as stack_models that facilitates rearranging the models into this flattened format.

Arguments

• models, data, logAge, metallicities are as in fit_sfh.
• Nsteps is the number of Monte Carlo samples you want to draw.

Keyword Arguments

• ϵ is the HMC step size. Convergence of the HMC samples is checked after sampling and if a convergence warning is issued, you should decrease this value.
• reporter is a valid reporter type from DynamicHMC.jl, either NoProgressReport, ProgressMeterReport for a basic progress meter, or LogProgressReport for more detailed reporting.
• show_convergence if true, will send sample convergence statistics to the default display.
• rng is a Random.AbstractRNG sampler instance that will be used when generating the random samples.

Returns

A NamedTuple with two elements:

• posterior_matrix is a Matrix with dimensions (npar, Nsteps) where npar is the number of fitting variables in the problem and is npar = length(bfgs_result.mle.μ). Each column is one independent sample.
• tree_statistics contains convergence statistics that can be viewed with DynamicHMC.Diagnostics.summarize_tree_statistics.

See also

• tsample_sfh for multi-threaded version.

tsample_sfh(bfgs_result::CompositeBFGSResult, 
            models::AbstractMatrix{<:Number},
            data::AbstractVector{<:Number},
            logAge::AbstractVector{<:Number},
            metallicities::AbstractVector{<:Number},
            Nsteps::Integer;
            ϵ::Real = 0.05, # HMC step size
            show_convergence::Bool=true,
            show_progress::Bool=true,
            rng::AbstractRNG=default_rng(),
            chain_length::Integer=100)
tsample_sfh(bfgs_result::CompositeBFGSResult, 
            models::AbstractVector{<:AbstractMatrix{<:Number}},
            data::AbstractMatrix{<:Number},
            logAge::AbstractVector{<:Number},
            metallicities::AbstractVector{<:Number},
            Nsteps::Integer;
            kws...)

Multi-threaded version of sample_sfh; see that method's documentation for details.

Implementation

This method splits the requested number of samples Nsamples into a number of independent HMC chains, each of which has length chain_length. Initial positions for each chain are randomly drawn from the multivariate Gaussian approximation to the objective function stored in bfgs_result, approximating a warm start. Smaller values of chain_length achieve better load balancing while larger values of chain_length allow each chain more time to mix (see also [8]). The default value of 100 results in good mixing with 24 fitting variables and a well-scaled step length ϵ – higher dimensional problems should increase chain_length. The downside to large chain_length is poor load balancing across available threads resulting in longer runtimes.

Notes

• if show_progress is true, we will show a progress bar that updates when individual chains complete. Currently this is not terribly useful unless the total number of chains is much greater than the number of available threads.

calculate_coeffs(MH_model::AbstractMetallicityModel,
                 disp_model::AbstractDispersionModel,
                 mstars::AbstractVector{<:Number}, 
                 logAge::AbstractVector{<:Number},
                 metallicities::AbstractVector{<:Number})

Returns per-SSP stellar mass coefficients ($r_{j,k}$ in the derivation) using the provided metallicity model MH_model and metallicity dispersion model disp_model for the set of SSPs with logarithmic ages logAge and metallicities metallicities.

Examples

julia> n_logage, n_mh = 10, 20; # Number of unique logAges, MHs

julia> coeffs = calculate_coeffs(PowerLawMZR(1.0, -1.0),
                                 GaussianDispersion(0.2),
                                 rand(n_logage),
                                 repeat(range(7.0, 10.0; length=n_logage); inner=n_mh),
                                 repeat(range(-2.0, 0.0; length=n_mh); outer=n_logage));

julia> coeffs isa Vector{Float64}
true

julia> length(coeffs) == n_logage * n_mh
true

(cum_sfh, sfr, mean_MH, samples, mstar, mstar_present) = 
    cum_sfr_quantiles(result::Union{CompositeBFGSResult,BFGSResult},
                      logAge::AbstractVector{<:Number},
                      MH::AbstractVector{<:Number}, T_max::Number,
                      Nsamples::Integer, q; surviving_fractions=nothing, kws...)

Draws Nsamples independent star formation histories from the solution result and calculates quantiles q across the samples in each unique time bin (i.e., unique(logAge)) for the cumulative star formation histories, star formation rates, and mean metallicities. Also calculates quantiles of the total birth stellar mass and, optionally, the present-day stellar mass of the sampled histories, and returns the drawn samples.

Arguments

• result::Union{CompositeBFGSResult,BFGSResult} is a BFGS result object as returned by fit_sfh, for example, whose contents will be used to sample random, independent star formation histories. If this is a CompositeBFGSResult containing both the MLE and MAP solutions, then the MLE solution is used for the best-fit values (SFRs and metallicity parameters) and the MAP solution is used for the uncertainty estimate.
• logAge::AbstractVector is a vector giving the log10(age [yr]) of the stellar populations that were used to derive result. For the purposes of calculating star formation rates, these are assumed to be left-bin edges.
• MH::AbstractVector is a vector giving the metallicities of the stellar populations that were used to derive result.
• T_max::Number is the rightmost final bin edge for calculating star formation rates in units of Gyr. For example, you might have logAge=[6.6, 6.7, 6.8] in which case a final logAge of 6.9 would give equal bin widths. In this case you would set T_max = exp10(6.9) / 1e9 ≈ 0.0079 so that the width of the final bin for the star formation rate calculation has the same log10(Age [yr]) step as the other bins.
• Nsamples::Integer is the number of random, independent star formation histories to draw from result to use when calculating quantiles.
• q are the quantiles you wish to have calculated. Can be one number (i.e., 0.5 for median), or most iterable types (e.g., (0.16, 0.5, 0.84) for median and 1-σ range).

Keyword Arguments

• surviving_fractions is nothing or a vector with one entry per stellar population (the same length as logAge) giving the fraction of its birth stellar mass in stars that are still alive, e.g., from surviving_mass_fraction. If provided, quantiles of the present-day stellar mass are returned as mstar_present.
• kws... are passed to calculate_cum_sfr, see that method's documentation for more information.

Returns

• cum_sfh::Matrix has size (length(unique(logAge)), length(q)) and contains the normalized cumulative SFH. This is ~1 at the most recent time in logAge and decreases as logAge increases.
• sfr::Matrix has size (length(unique(logAge)), length(q)) and gives the star formation rate across each bin in unique(logAge).
• mean_MH::Matrix has size (length(unique(logAge)), length(q)) and gives the stellar-mass-weighted mean metallicity of the stellar population as a function of unique(logAge).
• samples::Matrix has size (length(unique(logAge)), Nsamples) and contains the unique samples drawn from result that were used to derive the quantiles.
• mstar::Vector has length length(q) and gives the quantiles of the total birth stellar mass of the sampled star formation histories, in the units of the stellar mass coefficients of result. These are quantiles of the total mass of each sampled history, which differ from the sums of the quantiles of the individual time bins because the time bins are correlated.
• mstar_present is nothing if surviving_fractions is nothing, and otherwise has length length(q) and gives the quantiles of the present-day stellar mass (the sum over stellar populations of their stellar mass coefficients times surviving_fractions) of the same sampled star formation histories.

Examples

julia> result = fit_sfh(...)
CompositeBFGSResult{...} ...

julia> q = (0.16, 0.5, 0.84) # quantiles we want

julia> result = cum_sfr_quantiles(result, logAge, MH, 13.7, 10_000, q);

julia> length(result) == 6
true

julia> all(size(result[i]) == (length(unique(logAge)), length(q)) for i in 1:3)
true

julia> size(result[4], 2) == 10_000

tau(τ, unique_logAge, max_logAge, cum_sfh [, lower, upper])

Returns the lookback time (in Gyr) at which τ percent of the total stellar mass of a galaxy formed. See tau_interp for the underlying interpolator function.

Arguments

• τ: Fraction of a galaxy's total stellar mass for which you wish to know when the galaxy had that much mass (e.g., 0.50 for τ_50).
• unique_logAge should contain the list of unique log10(age [yr]) for which the SFH solution was derived.
• max_logAge should be the maximum log10(age [yr]) that you wish to consider – this sets when the stellar mass of the galaxy was 0.
• cum_sfh is an array of the cumulative stellar mass of the galaxy corresponding to the lookback times in unique_logAge; this can be calculated with calculate_cum_sfr or cum_sfr_quantiles.

Optional Arguments

• lower is an array of the cumulative stellar mass of the galaxy at some lower confidence level (e.g., 16% for 1-σ); if provided, the function will also return the lower bound on t.
• upper is an array of the cumulative stellar mass of the galaxy at some higher confidence level (e.g., 84% for 1-σ); if provided, the function will also return the upper bound on t.

Returns

• t::Number the lookback time in Gyr when the galaxy's stellar mass was τ times its total birth stellar mass.

If lower and upper are provided, a tuple of three numbers is returned: (t_lower, t_mle, t_upper).

Examples

First we will verify that tau returns the correct lookback time when the requested τ value corresponds exactly to a value in cum_sfh:

julia> unique_logAge = [8.0, 8.5, 9.0, 9.5, 10.0];

julia> max_logAge = 10.13;

julia> cum_sfh = [1.0, 0.8, 0.5, 0.2, 0.1];

julia> StarFormationHistories.tau(0.5, unique_logAge, max_logAge, cum_sfh) # Exact result
1.0

Now we will show an example where τ does not correspond exactly to a value in cum_sfh, in which case the function will interpolate between the two nearest values in cum_sfh to return the lookback time corresponding to the requested τ. The interpolation is performed in linear age space.

julia> cum_sfh = [1.0, 0.8, 0.4, 0.2, 0.1]; # Interpolation between log(age) = 8.5 and 9.0

julia> StarFormationHistories.tau(0.5, unique_logAge, max_logAge, cum_sfh) ≈ 0.8290569415042095
true

We can also pass multiple τ values at once:

julia> StarFormationHistories.tau([0.5, 0.75], unique_logAge, max_logAge, cum_sfh) ≈ [0.8290569415042095, 0.40169929526473325]
true

This will also work with reversed cum_sfh and unique_logAge:

julia> StarFormationHistories.tau(0.5, reverse(unique_logAge), max_logAge, reverse(cum_sfh)) ≈ 0.8290569415042095
true

Now we can pass in lower and upper estimates on the cumulative stellar mass to get confidence intervals on τ; the result is returned as a size (length(τ), 3) matrix where the first column is the lower bound on τ, the second column is the best estimate of τ, and the third column is the upper bound on τ:

julia> upper = [1.0, 0.85, 0.6, 0.3, 0.15];

julia> lower = [1.0, 0.75, 0.3, 0.1, 0.05];

julia> StarFormationHistories.tau(0.5, unique_logAge, max_logAge, cum_sfh, lower, upper) ≈ [0.6961012293408169  0.8290569415042095  1.7207592200561264]
true

And calling with a vector of τ values also works in this case:

julia> StarFormationHistories.tau([0.5, 0.75], unique_logAge, max_logAge, cum_sfh, lower, upper) ≈ [0.6961012293408169  0.8290569415042095  1.7207592200561264; 0.31622776601683794  0.40169929526473325  0.5897366596101027]
true

tau(result, τ, logAge, MH, max_logAge; Nsamples=10_000, q=(0.16, 0.5, 0.84), kws...)

Returns the lookback time (in Gyr) at which τ percent of the total stellar mass of a galaxy formed with upper and lower uncertainty estimates. Uses the solution result to draw samples of the cumulative star formation history and then calculates quantiles across those samples.

Arguments

• result::Union{CompositeBFGSResult, BFGSResult} is a BFGS result object as returned by fit_sfh, for example, whose contents will be used to sample random, independent star formation histories via cum_sfr_quantiles.
• τ: Fraction of a galaxy's total stellar mass for which you wish to know when the galaxy had that much mass (e.g., 0.50 for τ_50). Can also be a vector of multiple τ values (e.g., [0.5, 0.75, 0.9]).
• logAge::AbstractVector is a vector giving the log10(age [yr]) of the stellar populations that were used to derive result. For the purposes of calculating star formation rates, these are assumed to be left-bin edges.
• MH::AbstractVector is a vector giving the metallicities of the stellar populations that were used to derive result.
• max_logAge should be the maximum log10(age [yr]) that you wish to consider – this sets when the stellar mass of the galaxy was 0.

Keyword Arguments

• Nsamples::Integer is the number of random, independent star formation histories to draw from result to use when calculating quantiles.
• q is passed to cum_sfr_quantiles to calculate quantiles on the cumulative star formation history samples. This must be a length 3 tuple of quantiles (e.g., (0.16, 0.5, 0.84)) to get lower, best, and upper estimates on τ. It is not recommended that you change this.
• kws... are passed to calculate_cum_sfr, see that method's documentation for more information.

Returns

• t::Matrix, a size (length(τ), 3) matrix where the first column is the lower bound on τ, the second column is the best estimate of τ, and the third column is the upper bound on τ.

tau_interp(unique_logAge, max_logAge, cum_sfh)

Returns an interpolator for the lookback time (in Gyr) at which a given fraction of the total stellar mass of a galaxy formed. This is a helper function for tau.

AbstractBackground is the abstract supertype of the background models accepted by fit_sfh, NoBackground and HessBackground.

NoBackground()

Background model for fits without a background, so that every observed star must be explained by the stellar population templates. This is the default for fit_sfh.

HessBackground(shape::AbstractArray{<:Real}, b::Union{Real, Nothing}=nothing; floor::Real=0.05, free::Bool=true)

Background model for fit_sfh in which the expected number of stars in Hess diagram bin $i$ is

\[m_i = \sum_{j,k} r_{j,k} \, c_{i,j,k} + b \, u_i\]

where $u_i$ is the background shape and $b$ is the expected number of background stars in the fitted bins. Background stars are contaminants such as unresolved background galaxies or foreground stars, or stars the templates do not model. Without a background, observed stars in bins where the templates are (nearly) zero distort the fit, biasing the stellar mass high and the metallicity dispersion broad.

shape is the background Hess diagram over the same bins as the data passed to fit_sfh, e.g., the Hess diagram of a nearby control field, or ones(length(data)) for a uniform background. It is normalized to sum to 1, after which the fraction floor of the background is spread uniformly over the bins, so that stars in bins where both shape and the templates are zero can be explained by the background. b is the initial value of the expected number of background stars (a number of stars, not a fraction), which is fit if free is true and fixed otherwise. If b is nothing, fit_sfh starts from 1% of the observed stars. Like the stellar mass coefficients $R_j$, $b$ is fit in logarithmic coordinates, so it must be positive to be fit, and starting far from its best-fit value (e.g., at a single star when there are thousands of contaminants) slows or destabilizes the fit. A fixed b must be given, e.g., the number of stars in a control field scaled by the ratio of the areas.

A uniform background removes most of the bias from contamination spread over the Hess diagram, and costs nothing when there is no contamination, since $b$ then goes to zero. Contamination concentrated in part of the Hess diagram, such as faint background galaxies, should be modeled with its own shape; a uniform background explains only part of it.

After fitting, the best-fit background model is stored in the bg_model field of the BFGSResult, and its expected Hess diagram is bg_model.b .* bg_model.shape. $b$ is the last element of the fitted parameters μ.