Why is simulation used for probit choice probabilities? As explained in this answer, the probability of a certain outcome in a Logit model can be written as $$
P=\int_{\varepsilon=-\beta'x}^{\infty} f(\varepsilon)d\varepsilon\\
= 1- F(-\beta'x) = 1-\dfrac{1}{1+\exp(\beta'x)} = \dfrac{\exp(\beta'x)}{1+\exp(\beta'x)}
$$
I think I understand the reasoning behind this; the logit model assumes that the error term is distributed logistically, so the CDF of the logistic distribution is used. The probit model assumes that the error term follows a normal distribution giving
$$
 \Pr(Y=1 \mid X) = \Phi(X'\beta) =\int_{\varepsilon=-\beta'x}^{\infty} f(\varepsilon)d\varepsilon = \Phi(x_{i}'\beta) = \int_{-\infty}^{x_{i}'\beta} \dfrac{1}{\sqrt{2\pi}}e^{-\dfrac{1}{2}t^{2}}dt
$$
where $\Phi(x)$ is the CDF of the standard normal distribution. Now the referred answer, along with other sources, say that because these probabilities aren't closed-form, unlike those of the logit model, probit models require simulation.
But why use something complicated as simulation when there are a variety of methods that very accurately approximate $\Phi(x)$? It would seem to me that this would both be less complicated and more accurate.
 A: For a simple probit, software I am familiar with (Stata, R, SPSS) does not use simulation to compute that integral. 
Simulation is sometimes used in more complicated variations, like a multinomial probit (this is what the answer you linked to refered to) or when random effects are included in the model. There the likelihood function contains an additional integral over the additional (group level) error term(s). For a list of examples see this special issue of the Stata Journal.
A: In both cases the model is
$$Y\sim Be(p)$$
in which Y is the dependent variable. Here we model the effect of independent variables on Y through a link function. Logit and probit functions are used as the link function for logisic and probit regression models respectively
 $$logit(p)=X'\beta~~~~or~~~~~\Phi(p)=X'\beta$$
So it's easy to see if you want to simulate Y based on an observation of X, you have to compute p based on the $$X'\beta$$ In logistic regression it is straight forward as you have mentioned. But in probit regression model you have to compute $$p=\Phi^{-1}(X'\beta)$$ to be able to simulate Y. So it's up to you which method you like to use to compute this quantity. You can use numerical methods or interpolating from a table or any other method you know.
