Probit model question (regression) I'm reading a thesis and I need your help to understand the equation below.  
$$\Pr(\text{failure}=1 \mid X_1,X_3,X_3,X_4)=\int_{-\infty}^z \varphi(k) \, dk\tag{1}$$
$\varphi(k)$ is the standard normal density and $z=\beta_1 X_1+\beta_2 X_2+\beta_3 X_3+\beta_4 X_4$
We estimate regression (1) through probit.   
It means,
first we estimate probability of failure ($z$, in here) by probit test,
and then do the normal regression to estimate $\beta_1,\beta_2, \beta_3,\beta_4$?
(by using $z$ from first stage and we already have $X_1,X_2,X_3,X_4$ values)
We have to regress twice in this model?
 A: In a probit model the data would typically consist of a list of cases, and for each case one would have observed the values of $X_1,X_2,X_3,X_4$ (the number of such predictors may vary) and whether there was a success or a failure.  Let $n$ be the number of cases.  One observes $(I_i,x_{1,i},\ldots,x_{4,i})$ for $i=1,\ldots,n$, and the value of $I_i$ would be "success" or "failure", and often this would be represented by $1$ or $0$, the number of successes in the $i$th case.  But below I'll do it the other way, for sake of consistency with your notation, i.e.
\begin{align}
0 & = \text{success}, \\ 1 & = \text{failure}.
\end{align}
One estimates $\beta_1,\beta_2,\beta_3,\beta_4$ by maximum likelihood.  That means the estimates $\widehat\beta_1,\widehat\beta_2,\widehat\beta_3,\widehat\beta_4$ would be the values that maximize the probability of the data that were actually observed.  To be more explicit: The likelihood function is
$$
L(\beta_1,\beta_2,\beta_3,\beta_4) = \left(\prod_{i\,:\,I_i=1} \int_{-\infty}^{\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}} \varphi(k)\,dk \right)\left(\prod_{i\,:\,I_i=0} \left(1 - \int_{-\infty}^{\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}} \varphi(k)\,dk\right) \right).
$$
Each integral is from $-\infty$ to $\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}$.
The estimates $\widehat\beta_1,\widehat\beta_2,\widehat\beta_3,\widehat\beta_4$ are the values that maximize this.  These would be found by numerical methods; I don't think there's a closed form.
I think typically there would be an additional parameter: one would have $\alpha+\beta_1 x_{1,i}+\cdots+\beta_4 x_{4,i}$ as the upper bound of integration, and the likelihood function would be $L(\alpha,\beta_1,\beta_2,\beta_3,\beta_4)$.
So: No, you don't regress twice.
