I am taking first steps with Julia, trying to use Optim.jl to estimate a probit. I have adapted some code I [found] :

```
using DataFrames
using Distributions
using Optim
# -------------- Normal distribution functions prevent 0.0 and 1.0 ----------------------
function normcdf(x::Vector{Float64})
out = Distributions.cdf(Distributions.Normal(), x)
out + (out .== 0.0)*eps(1.0) - (out .== 1.0)*eps(1.0)
end
function normpdf(x::Vector{Float64})
Distributions.pdf(Distributions.Normal(), x)
end
# -------------- Likelihood and gradient ----------------------
function λ(β::Vector{Float64})
q = 2d-1
q .* normpdf(q .* X*β) ./ normcdf(q.*X*β)
end
function LL(β::Vector{Float64})
- sum( log.( normcdf( (2d-1) .* X*β) ) )
end
function g!(storage, β)
storage[:] = -sum( λ(β).*X, 1 )
end
# ------------------------ Probit DGP ------------------------
N = 10^5
k = 4
u = rand(Normal(),N)
X = [ones(N) reshape( rand(Uniform(-1,1), N*(k-1)), (N,k-1) )]
βTrue = [2.0, 3, 5, 10]
y = X*βTrue + u
d = 1.0 * (y .>= 0 )
```

I then run this, first without supplying the gradient function:

```
res1 = Optim.optimize(LL, zeros(k), method=BFGS())
```

which works fine but when I supply the gradient function:

```
res2 = Optim.optimize(LL, g!, zeros(k), method=BFGS())
```

it reports convergence despite not iterating and the results being clearly incorrect. I presume I have made a simple mistake but I can’t see it.