We consider stochastic gradient estimation using only black-box function evaluations, where the function argument lies within a probability simplex. This problem is motivated from gradient-descent optimization procedures in multiple applications in distributionally robust analysis and inverse model calibration involving decision variables that are probability distributions. We are especially interested in obtaining gradient estimators where one or few sample observations or simulation runs apply simultaneously to all directions. Conventional zeroth-order gradient schemes such as simultaneous perturbation face challenges as the required moment conditions that allow the "canceling" of higher-order biases cannot be satisfied without violating the simplex constraints. We investigate a new set of required conditions on the random perturbation generator, which leads us to a class of implementable gradient estimators using Dirichlet mixtures. We study the statistical properties of these estimators and their utility in constrained stochastic approximation, including both Frank-Wolfe and mirror descent update schemes. We demonstrate the effectiveness of our procedures and compare with benchmarks via several numerical examples.