Improve the residual calculation for regmodsm

[AHsu98]

In the get_residual_computation_function we're defining lambda functions for all the residuals individually.

There's a todo comment that we want to vectorize it, that shouldn't be too hard, but depends on how row is passed in later, may require refactoring downstream usage.

We also don't have a description for the string 'sigma' parameter in the docstring.

We're taking a complex step for the derivative of tobit, this should be completely avoidable directly differentiating the gaussian CDF.

I'd propose we abstract to the calculus, and call subroutines for the different likelihoods and links, but I need to better understand the downstream usage better before being able to implement a vectorized version.