mata:
void plleval3(real scalar todo, real vector b,     ///
              real vector y,    real matrix X,     ///
	      val, grad, hess)
{
    real vector  xb

    xb = X*b'
   val = -exp(xb) + y:*xb - lnfactorial(y)
}

y  = st_data(., "accidents")
X  = st_data(., "cvalue ibn.kids traffic")
X  = X,J(rows(X), 1, 1)

C  = e(5, 7)
c  = 0
Ct = C,c

S  = optimize_init()
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_evaluator(S, &plleval3())
optimize_init_evaluatortype(S, "gf0")
optimize_init_params(S, J(1, 7, .01))
optimize_init_constraints(S, Ct)

bh = optimize(S)
optimize_result_params(S)
end