GitHub

Metallicity Dispersion Models

In order to reproduce the broadness of features in observed color-magnitude diagrams it is common to introduce some dispersion in metallicity for stars formed in each time bin, where each time bin is associated with a mean metallicity $\mu_j$. We implement a generic hierarchical model to accomodate this functionality.

The hierarchical model forms the per-SSP weights $r_{j,k}$, which are indexed by population age $j$ and metallicity $k$, as a function of a linear coefficient $R_j$ which describes the stellar mass formed in the time bin, and a relative weight $A_{j,k}$ which depends on the mean metallicity $\mu_j$ and the metallicity of the SSP under consideration. In the case of a Gaussian metallicity dispersion at fixed age, which is often used in practice, we can write

\[A_{j,k} = \text{exp} \left( - \frac{1}{2} \left( \frac{ [\text{M}/\text{H}]_k - \mu_j}{\sigma} \right)^2 \right)\]

If we take $R_j$ to be the total stellar mass formed in the time bin, then it is clear that we require the sum over the relative weights for the SSPs of age $j$ to equal one, i.e., $\sum_k r_{j,k} = 1$. We therefore require that the relative weight on each SSP template of age $j$ be normalized by the sum $\sum_k A_{j,k}$, so that the relative weights are

\[r_{j,k} = R_j \, \frac{A_{j,k}}{\sum_k A_{j,k}}\]

We provide a generic interface for describing the analytic form of the $A_{j,k}$ so that it is easy to define new dispersion models that will integrate with our fitting routines. Built-in, ready to use models are described below, and the API for defining new models is described in the API section.

Built-In Models

StarFormationHistories.GaussianDispersion — Type
GaussianDispersion(σ::Real, free::NTuple{1, Bool} = (true,)) <: AbstractDispersionModel

Dispersion model for a Gaussian (i.e., Normal) spread in metallicities with standard deviation σ (which must be greater than 0) at fixed age. The relative weights for this model are given by $A_{j,k} = \exp(-((x_k - μ_j)/σ)^2/2).$ The σ can be fit during optimizations if free == (true,) or fixed if free == (false,).

Examples

julia> GaussianDispersion(0.2) isa GaussianDispersion{Float64}
true

julia> import Test

julia> Test.@test_throws(ArgumentError, GaussianDispersion(-0.2)) isa Test.Pass
true

julia> nparams(GaussianDispersion(0.2)) == 1
true

julia> GaussianDispersion(0.2)(1.0, 1.2) ≈ exp(-0.5)
true

julia> all(values(gradient(GaussianDispersion(0.2), 1.0, 1.2)) .≈
                (3.0326532985631656, -3.0326532985631656))
true

julia> update_params(GaussianDispersion(0.2), 0.3) == GaussianDispersion(0.3)
true

julia> transforms(GaussianDispersion(0.2)) == (1,)
true

julia> free_params(GaussianDispersion(0.2, (false,))) == (false,)
true
source

Metallicity Dispersion API

Below we describe the API that must be followed in order to implement new types for describing the $A_{j,k}$, such that they will work with our provided fitting and sampling methods.

StarFormationHistories.AbstractDispersionModel — Type

Abstract type for all models of metallicity dispersion at fixed time $t_j$, for which the mean metallicity is $\mu_j$. Concrete subtypes T <: AbstractDispersionModel should implement the following API:

  • (model::T)(x::Real, μ::Real) should be defined so that the struct is callable with a metallicity x and a mean metallicity μ, returning the relative weight for the metallicity x given the dispersion model. This is $A_{j,k}$ for $\mu = \mu_j$ in the derivations presented in the documentation.
  • nparams(model::T) should return the number of fittable parameters in the model.
  • fittable_params(model::T) should return the values of the fittable parameters in the model.
  • gradient(model::T, x::Real, μ::Real) should return a tuple that contains the partial derivative of the $A_{j,k}$ with respect to each fittable model parameter, plus the partial derivative with respect to μ as the final element.
  • update_params(model::T, newparams) should return a new instance of T with the fittable parameters contained in newparams (which is typically a vector or tuple) and non-fittable parameters inherited from the provided model.
  • transforms(model::T) should return a tuple of length nparams(model) which indicates how the fittable variables should be transformed for optimization, if at all. Elements should be 1 for parameters that are constrained to always be positive, 0 for parameters that can be positive or negative, and -1 for parameters that are constrained to always be negative.
  • free_params(model::T) should return an NTuple{nparams(model), Bool} that is true for fittable parameters that you want to optimize and false for fittable parameters that you want to stay fixed during optimization.
source
StarFormationHistories.gradient — Method
gradient(model::AbstractDispersionModel{T}, x::Real, μ::Real)::NTuple{nparams(model)+1, T}

Returns a tuple containing the partial derivative of the model with respect to all fittable parameters, plus the partial derivative with respect to the mean metallicity μ as the final element. These partial derivatives are evaluated at metallicity x where the model has expectation value μ.

source
StarFormationHistories.update_params — Method
update_params(model::T, newparams)::T where {T <: AbstractDispersionModel}

Returns a new instance of the model type T with the fittable parameters contained in newparams (which is typically a vector or tuple), with non-fittable parameters inherited from the provided model.

source
StarFormationHistories.transforms — Method
transforms(model::AbstractDispersionModel)::NTuple{nparams(model), Int}

Returns a tuple of length nparams(model) which indicates how the fittable variables should be transformed for optimization, if at all. Elements should be 1 for parameters that are constrained to always be positive, 0 for parameters that can be positive or negative, and -1 for parameters that are constrained to always be negative.

source
StarFormationHistories.free_params — Method
free_params(model::AbstractDispersionModel)::NTuple{nparams(model), Bool}

Returns an tuple of length nparams(model) that is true for fittable parameters that you want to optimize and false for fittable parameters that you want to stay fixed during optimization.

source