function (DF) for a single star is denoted $f_{\bm{W}}(\bm{w})$ [for an extensive discussion, see @binney_galactic_2008].[^2] We assume that the DF is an element of the model $\{f(\bm{w};\bm{a})\}_{\bm{a} \in \bm{A}}$, where the parameters are $D$-dimensional real vectors, the elements of which represent the total galactic mass, galactic scale length, velocity anisotropy, etc. From the DF we may calculate the observable properties of the system. For a dSph these observables are typically the projected stellar positions, (represented by the random variables $X$ and $Y$) and the line-of-sight velocity (represented by the random variable $V_z$). The DF for these observables is given by the marginalization of the phase-space DF: $$\begin{aligned}
    f_{\bm{W}_{\mathrm{p}}}(\bm{w}_{\mathrm{p}}; \bm{a}) = \int_{\mathbf{R}^3}
    f_{\bm{W}}(\bm{w}; \bm{a}) {\operatorname{d}}z {\operatorname{d}}v_{x} {\operatorname{d}}v_{y}.\end{aligned}$$ We can reasonably assume that the states of stars are independent and identically distributed [@binney_galactic_2008] so that the joint marginalized PDF for $N$ stars is $$\begin{aligned}
    f_{(\bm{W}_{\mathrm{p}, 1}, \ldots, \bm{W}_{\mathrm{p},
        N})}(\bm{w}_{\mathrm{p}, 1}, \ldots, \bm{w}_{\mathrm{p}, N}; \bm{a}) =
    \prod_{i = 1}^{N} f_{\bm{W}_{\mathrm{p}}}(\bm{w}_{\mathrm{p}, i}; \bm{a}).\end{aligned}$$ By comparing these observables with a data set, we can then optimize the parameters of the phase-space DF using the maximum likelihood method.

By definition, the likelihood of model parameter $\bm{a}$ is $$\begin{aligned}
    L(\bm{a}; \bm{w}_{\mathrm{p}, 1}, \ldots, \bm{w}_{\mathrm{p}, N}) =
    f_{(\bm{W}_{\mathrm{p}, 