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?