  We present a method for fitting orbit-superposition models to the kinematics
of discrete stellar systems when the available stellar sample has been filtered
by a known selection function. The fitting method can be applied to any model
in which the distribution function is represented as a linear superposition of
basis elements with unknown weights. As an example, we apply it to Fritz et
al.'s kinematics of the innermost regions of the Milky Way's nuclear stellar
cluster. Assuming spherical symmetry, our models fit a black hole of mass
$M_\bullet=(3.76\pm0.22)\times10^6\,M_\odot$, surrounded by an extended mass
$M_\star=(6.57\pm0.54)\times10^6\,M_\odot$ within 4 pc. Within 1 pc the
best-fitting mass models have an approximate power-law density cusp
$\rho\propto r^{-\gamma}$ with $\gamma=1.3\pm0.3$. We carry out an extensive
investigation of how our modelling assumptions might bias these estimates:
$M_\bullet$ is the most robust parameter and $\gamma$ the least. Internally the
best-fitting models have broadly isotropic orbit distributions, apart from a
bias towards circular orbits between 0.1 and 0.3 parsec.
