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.