mata: void MYNLProbit(real scalar todo, real vector b, /// real vector y, real matrix X, /// val, grad, hess) { real vector r, f, xb real matrix df xb = X*b' f = normal(xb) r = y - f val = -(r:^2) df = normalden(xb):*X if (todo>=1) { grad = r:*df } if (todo==2) { hess = -1*quadcross(df, df) } } y = st_data(., "hadaccident") X = st_data(., "cvalue tickets") n = rows(y) X = X, J(n, 1, 1) p = cols(X) S = optimize_init() optimize_init_argument(S, 1, y) optimize_init_argument(S, 2, X) optimize_init_evaluator(S, &MYNLProbit()) optimize_init_conv_nrtol(S, 1e-11) optimize_init_params(S, J(1, 3, .01)) optimize_init_evaluatortype(S, "gf2") bh = optimize(S) M = invsym(-1*optimize_result_Hessian(S)) sb = (-1/(n-p))*optimize_result_value(S) V = sb*M "Point estimates" bh "VCE" V end