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 Ways nuclear stellar cluster. Assuming spherical symmetry, our models fit a black hole of mass $M_bullet=(3.76pm0.22)times10^6,M_odot$, surrounded by an extended mass $M_star=(6.57pm0.54)times10^6,M_odot$ within 4 pc. Within 1 pc the best-fitting mass models have an approximate power-law density cusp $rhopropto r^{-gamma}$ with $gamma=1.3pm0.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.