Probit model with Gaussian noise

Question

Assume we have the following model setup $$\Phi^{-1}(D)=\alpha+\beta X+\epsilon$$ where $\epsilon\sim N(0,\sigma^{2})$ and $D_{i}=\{0,1\}$. This implies that $$\text{Pr}(D_{i}=1\,|\,X,\epsilon)=\Phi(\alpha+\beta X+\epsilon)$$ but we need to integrate out the $\epsilon$ $$\begin{align} \text{Pr}(D_{i}=1\,|\,X)&=\int_{\epsilon}\text{Pr}(D_{i}=1\,|\,X,\epsilon)\,f_{\epsilon}(\epsilon)\,d\epsilon\\ &=\Phi\bigg(\frac{\alpha+\beta x}{\sqrt{1+\sigma^{2}}}\bigg) \end{align}$$

So if we were to estimate the parameters of this model $(\hat{\alpha},\hat{\beta},\hat{\sigma})$ we could go about it by changing the likelihood function of the simple probit model from $$\begin{align} L=\sum_{i=1}^{n}D_{i}\log(\Phi(\alpha+\beta X+\epsilon))+(1-D_{i})\log(\Phi(\alpha+\beta X+\epsilon)) \end{align}$$ to $$\begin{align} L^{*}=\sum_{i=1}^{n}D_{i}\log\bigg(\Phi\bigg(\frac{\alpha+\beta X}{\sqrt{1+\sigma^{2}}}\bigg)\bigg)+(1-D_{i})\log\bigg(1-\Phi\bigg(\frac{\alpha+\beta X}{\sqrt{1+\sigma^{2}}}\bigg)\bigg) \end{align}$$ However, I've noticed that performing reliable optimisation of this likelihood is difficult. Given the toy example

n = 10000
a = -2
b = 0.01
x = runif(n, min = 1, max = 5) + rnorm(n, 0, 0.15)
p = pnorm(a + b*x)
d = rbinom(n, size = 1, prob = p)
y = tibble::as_tibble(data.frame(x, p, d))

and the likelihood defined as

fn = function(par, x, d) {
  return(-sum(
    d*log(pmax(10^-23, pnorm((par[1] + par[2]*x)/sqrt(1 + par[3]^2)))) +
      (1-d)*log(pmax(10^-23, 1 - pnorm((par[1] + par[2]*x)/sqrt(1 + par[3]^2))))
  ))
}

and using quasi-Newton methods

optim(par = c(0, 0, 0.5),
      fn = fn,
      x = y$x,
      d = y$d,
      method = "L-BFGS-B",
      lower = c(-Inf, -Inf, 0),
      upper = c(Inf, Inf, Inf),
      hessian = TRUE)

typically doesn't behave very well. In fact, the $\hat{\sigma}$ usually just converges to a point near the starting value.

Are there any obvious changes (choice of algorithm, approximations to the likelihood function, better choice of starting values) that can be made to make the estimation of $(\hat{\alpha},\hat{\beta},\hat{\sigma})$ more reliable?

Answer 1

The model is not identified, meaning there is no unique solution to the optimization problem. There are infinite values of the parameters that will yield the same likelihood. For example, $\alpha = .5$, $\beta = 1$, and $\sigma = 2$ will yield the exact same likelihood as $\alpha = 1$, $\beta = 2$, and $\sigma = \sqrt{19}$.

More generally, consider the maximum of the likelihood, $L^*$, which is found when $(\alpha, \beta, \sigma) = (\alpha^*, \beta^*, \sigma^*)$. For any $k$, $$\left(k\alpha^*, k\beta^*, \sqrt{k^2+1+(k\sigma^*)^2}\right)$$ will yield the exact same likelihood. Therefore, there is no unique value of the parameters that maximizes the likelihood. This is why the optimization is unstable; any specific solution it arrives at will be due purely to numerical instability.

Note that in typical probit regression, we assume $\sigma=0$, i.e., that there is no latent variable $\epsilon$ that is unaccounted for. This is different from the latent variable formulation of probit regression, where we assume $$Y^*=X\beta+\epsilon$$ where $\epsilon \sim N(0, 1)$, and $$P(D=1|X) = P(Y^*>0|X)=P(X\beta + \epsilon > 0)$$ which implies $P(D=1|X) = \Phi(X\beta)$.