We describe a maximum-likelihood method for determining the mass distribution in spherical stellar systems from the radial velocities of a population of discrete test particles .
The method assumes a parametric form for the mass distribution and a non-parametric two-integral distribution function .
We apply the method to a sample of 161 globular clusters in M87 .
We find that the mass within 32 kpc is ( 2.4 \pm 0.6 ) \times 10 ^ { 12 } M { { } _ { \odot } } , and the exponent of the density profile \rho \propto r ^ { - \alpha } in the range 10 - 100 kpc is \alpha = 1.6 \pm 0.4 .
The energy distribution suggests a few kinematically distinct groups of globular clusters .
The anisotropy of the globular-cluster velocity distribution can not be determined reliably with the present data .
Models fitted to an NFW potential yield similar mass estimates but can not constrain the concentration radius r _ { c } in the range 10 - 500 kpc .