We investigate the dark matter density profile of the massive elliptical galaxy, NGC 1407, by constructing spherically symmetric Jeans models of its field star and globular cluster systems. Two major challenges in such models are the degeneracy between the stellar mass and the dark matter halo profiles, and the degeneracy between the orbital anisotropy of the tracer population and the total mass causing the observed motions. We address the first issue by using new measurements of the mass-to-light ratio profile from stellar population constraints that include a radially varying initial mass function. To mitigate the mass-anisotropy degeneracy, we make use of multiple kinematic tracers, including two subpopulations of globular clusters in addition to the galaxys field stars. We create a hierarchical Bayesian model that addresses several often-neglected systematic uncertainties, such as the statistical weight given to various data sets and the adopted distance. After sampling the posterior probability distribution with a Markov chain Monte Carlo method, we find evidence for a central cusp with a log slope of γ = 1.0-0.4 +0.2 (stat) -0.5 0.3 (sys), with the quantified systematic uncertainty dominated by choice of anisotropy profile. This is lower than expected for dark matter halos that have undergone adiabatic contraction, supporting inferences from gravitational lensing that some process has suppressed the steepening of halos in massive galaxies. We also confirm radially biased orbits for the metal-rich globular clusters and tangentially biased orbits for the metal-poor globular clusters, which remains a puzzling finding for an accretion-dominated halo.