Major update. Implemented hessian matrix for family = ‘binomial’. Much
faster than previous version where numerical derivative was used.