We apply a new simulation method that solves the multidimensional probability integrals that arise in maximum likelihood estimation of a broad class of limited dependent variable models. The simulation method has four key features: the simulated choice probabilities are unbiased; they are a continuous and differentiable function of the parameters of the model; they are bounded between 0 and 1; and their computation takes an effort that is nearly linear in the dimension of the probability integral, independent of the magnitudes of the true probabilities. We also show that the new simulation method produces probability estimates with substantially smaller variance than those generated by acceptance-rejection methods or by Stern’s (1987) method. The simulated probabilities can therefore be used to revive the Lerman and Manski (1981) procedure of approximating the likelihood function using simulated choice probabilities by overcoming its computational disadvantages.