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.

We need a one-line description at the top of classes, methods, and functions for the docstring rendering. https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt#sections