BUG+ENH: stats:  Improve the 'fit' method of norm_gen, gamma_gen and beta_gen.

This should close gh-2514.

The specific changes are summarized in the commit message.

An example of the improved performance:

```
In [5]: from scipy.stats import beta, rv_continuous

In [6]: x = beta.rvs(0.5, 1.5, 0, 4, size=500)

In [7]: beta.fit(x, floc=0, fscale=4)
Out[7]: (0.51490904775491819, 1.4783326683085065, 0, 4)

In [8]: rv_continuous.fit(beta, x, floc=0, fscale=4)
/home/warren/local_scipy/lib/python2.7/site-packages/scipy/stats/distributions.py:2442: RuntimeWarning: invalid value encountered in sqrt
  sk = 2*(b-a)*sqrt(a + b + 1) / (a + b + 2) / sqrt(a*b)
Out[8]: (0.51491189692021333, 1.4783061325955025, 0, 4)

In [9]: %timeit beta.fit(x, floc=0, fscale=4)
1000 loops, best of 3: 370 us per loop

In [10]: %timeit rv_continuous.fit(beta, x, floc=0, fscale=4)
100 loops, best of 3: 9.59 ms per loop

In [11]: 9.59/.37
Out[11]: 25.91891891891892
```
In this case, the new `fit` method is almost 26 times faster than `rv_continuous.fit`.