Re: Trying to understand Distributions

On Mon, Jul 12, 2021 at 4:06 PM Keith Sloan <[email protected]> wrote:

> Okay I changed the plot histogram as advised by Robert Kern ax2 as opposed
> to ax1
>
> but how do I get the correct gamma plot as I am not understanding
> np.linespace other the first parameter is a minimum and second maximum
>
> from astropy.table import Table, join
> import numpy as np
> import matplotlib.pyplot as plt
> from matplotlib.pyplot import plot
>
> RawMassEClassEmeasure =
> Table.read('../../GAMA_Data/REMassEClassEmeasure.fits')
> #print(RawMassEClassEmeasure.colnames)
> # CLEAN DATA
> #REMassEClassEmeasure =
> RawMassEClassEmeasure[RawMassEClassEmeasure['CountInCyl']> -500]
> RErange = RawMassEClassEmeasure[RawMassEClassEmeasure['CountInCyl']> -500]
> RErange1 = RErange[RErange['SurfaceDensity']< 50]
>
> binCount = 30
> alphaVal = .3
>
> ##### uminusr
> fig = plt.figure(figsize=(12, 6), dpi=200)
> fig.suptitle('Plot - Histogram Red Galaxies for Elliptical Galaxies')
> #fig.legend(loc="upper right")
> #import scipy.stats as stats
> from scipy import stats
> xfield = 'uminusr'
> counts, bins = np.histogram(RErange1[xfield].data,bins=binCount)
> print(counts)
> ag, bg, cg =stats.gamma.fit(counts)
> print(ag, bg, cg)
>

You don't use fit() on histogram counts. Use it on the data itself.

You probably want to fix the `loc` parameter when you do the fitting. It is
not typical to allow a zero-shifted gamma distribution (though I am told
that astronomy is a place where that does happen). We can use floc=0 to fix
that value and we can ignore it when it comes out of the output.

gam_shape, _, gam_scale = stats.gamma.fit(RErange1[xfield].data, floc=0.0)

Here's a quick demonstration with 100 points that I randomly draw from a
gamma distribution:

|3> data = stats.gamma.rvs(3, scale=2.0, size=100)


|4> data
array([ 3.07534256,  4.57813927, 13.8711741 ,  6.57534761,  5.28268719,
 7.7073074 ,  2.95929127,  7.58245554,  6.10107878,  3.52840919,
 2.91885749,  4.05885291,  5.03287215, 10.76342819,
        8.6441407 ,  8.64011327,  5.461707  ,  6.15034854,  2.84274024,
 8.16282873,  5.111948  , 11.0312346 ,  0.50150262,  4.31496838,
 5.78209971,  9.0395057 ,  3.97753423,  2.81673685,
        1.71633995,  9.24401038,  7.50455745,  3.90409068,  4.56853541,
16.55915147,  2.33004283,  2.89232249, 10.0872901 ,  3.67701972,
 9.15203612,  3.33052913,  5.38745436,  3.95826717,
        3.69810183,  5.52858701,  6.04013903,  1.41587988,  4.53861925,
 7.65106151,  4.05856947,  6.02678829,  6.08308182,  4.22620438,
 9.93515377,  9.45789334,  3.17088024, 10.59471748,
       12.19263839,  8.12720946,  8.05421344,  5.31204559,  4.65193752,
 6.83243734,  0.85212066, 16.05784841,  1.97212606,  6.28613507,
 8.00692823,  4.0366943 ,  5.07313536,  9.36507581,
        7.28824511,  6.38240736,  3.4051997 ,  6.24772202,  3.11340298,
 9.00223885,  6.87387583,  4.75206877,  3.68490716,  6.71069272,
 2.43124401, 11.01701724,  5.33665806,  4.36158851,
        5.89422855,  6.60390376,  8.24682823,  4.60149744,  5.11966798,
 1.72012758,  3.85308752,  4.13263933,  6.71598009,  5.47673371,
 2.89124835,  4.88451718,  6.69019804,  2.82428072,
        6.06868802, 11.30380543])

|5> stats.gamma.fit(data, floc=0)
(3.6909888310977736, 0, 1.6193634011500304)

ax2 = fig.add_subplot(3, 1, 2)
> ax2.hist(RErange1[xfield].data, bins=binCount, density=True)
> x = np.linspace(stats.gamma.ppf(0.1, ag),stats.gamma.ppf(0.99, ag), 243)
>
#ax2.plot(x, stats.gamma.pdf(x, ag),'r-', lw=5, alpha=0.6, label='gamma
> pdf')
>
Now be sure to use all of the fitted parameters, not just the shape
parameter (we can ignore the loc= parameter since we fixed it to 0).

x0, x1 = stats.gamma.ppf([0.1, 0.99], gam_shape, scale=gam_scale)
x = np.linspace(x0, x1, 243)
# Though this might be better:
# x = np.linspace(0, RErange1[xfield].data.max(), 243)
ax2.plot(x, stats.gamma.pdf(x, gam_shape, scale=gam_scale), ...)

-- 
Robert Kern

_______________________________________________
SciPy-User mailing list
[email protected]
https://mail.python.org/mailman/listinfo/scipy-user