In case a score method is not provided, this function will be used.
It computes the score by numerical differentiation of the log-likelihood.
The likelihood model
Custom arguments to pass to numDeriv::grad.
Additional arguments (to pass into loglik)