ENH: Add gradient-based optimization for fitting GLM's

closes #1736 and #1958, see also #1788 

This is a PR to allow GLM's to be fit using scipy optimizers.  It also fills some other gaps and corrects a few inconsistencies in the GLM implementation.  Here are the main points:

* The method keyword argument to fit can be IRLS (as before), or the name of any scipy optimizer.

* Several deriv2 implementations that were missing in family.py were falling back to numerical differentiation, which leads to issues with run time and big inaccuracies in some cases.  I have added most of the missing analytic deriv2 functions.

* I added tests for most of the links and associated functions (inverses, derivatives, etc.).

* I resolved inconsistencies in some of the signatures for the negative binomial family.

* I added a likelihood function that takes params as an argument, renaming the likelihood that takes the fitted means as loglike_mu.

* I added tests comparing IRLS fits to Newton fits for 6 families.  The tests all run clean for me.

Thanks, but they still need to be 79 characters. You can make a small one-line description. Then expand if need be. [Example](https://github.com/numpy/numpy/blob/master/numpy/linalg/linalg.py#L296)

Not trying to be picky. I just have to fix these things later when the docs build complains.