ENH: stats: Some enhancements to the chisquare functions:

* Add the `axis` argument to stats.chisquare and stats.mstats.chisquare.
* Add the `ddof` argument to stats.mstats.chisquare.

Closes gh-2108 (trac #1583).

The changes should be backwards compatible.

There is a change in behavior that I believe is a bug in the existing implementation.  Consider

    >>> f_obs = np.array([12, 12, 24, 12])
    >>> f_exp = np.array([[16, 16, 16, 12], [12, 12, 24, 12], [15, 15, 15, 15]])

In the call `chisquare(f_obs, f_exp)`, `f_obs` and `f_exp` are broadcast (which in this case results in an array with shape (3, 4)), and then a reduction operation computes the chi-squared statstic.  In the 0.12.0 implementation, the default is `axis=0`, which means the reduction that computes the chi-squared statistic operates on the wrong axis.  That is, in scipy 0.12.0 we get:

    >>> chisquare(f_obs, f_exp)
    (array([ 1.6,  1.6,  9.4,  0.6]), array([ 0.65938982,  0.65938982,  0.02441934,  0.89643237]))

Note that the arrays have shape (4,).  That is not the intended result.

The code in the pull request handles this case correctly:

    >>> chisquare(f_obs, f_exp)
    (array([ 6. ,  0. ,  7.2]), array([ 0.11161023,  1.        ,  0.06578905]))