We present a numerical method which is used to calculate anomalous scaling exponents of structure functions in the Kraichnan passive scalar advection model [R. H. Kraichnan, Phys. Fluids 11, 945 (1968)]. This Monte Carlo method, which is applicable in any space dimension, is based on the Lagrangian path interpretation of passive scalar dynamics, and uses the recently discovered equivalence between scaling exponents of structure functions and relaxation rates in the stochastic shape dynamics of groups of Lagrangian particles. We calculate third and fourth order anomalous exponents for several dimensions, comparing with the predictions of perturbative calculations in large dimensions. We find that Kraichnan's closure appears to give results in close agreement with the numerics. The third order exponents are compatible with our own previous nonperturbative calculations.