  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\pm0.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\pm0.4$.The energy distribution suggests a few kinematically
distinct groups of globular clusters. The anisotropy of the globular-cluster
velocity distribution cannot be determined reliably with the present data.
Models fitted to an NFW potential yield similar mass estimates but cannot
constrain the concentration radius $r_c$ in the range 10-500 kpc.
