On Mon, Jul 27, 2009 at 02:33, Daniel J Farrell<daniel.farrell@imperial.ac.uk> wrote:
Dear list,
I am looking for some of the functionality provided by the GNU Scientific Library histogram module (http://www.gnu.org/software/gsl/manual/html_node/Histograms.html ).
In particular, a need to be able to create a probability density function from my histogram of data. This will allow inverse look-ups to be performed, i.e. for a random number (0-->1) find the associated probability (http://www.gnu.org/software/gsl/manual/html_node/The-histogram-probability-d... ). This allows the distribution and for samples to be returned weighted by the probability of the distribution -- which is a common task!
It looks like you want a CDF (or rather it's inverse, the PPF) rather than a PDF. Anyways, this is straightforward. Compute the histogram using normed=False. Find the cumulative sum and divide by the sum to get the (smoothed) empirical CDF. Prepend a 0.0 to this, and then this will align with the edges array that is also returned. Then you can use linear interpolation to do lookups. If you use the edges array as "X" and the empirical CDF as "Y", then this is a CDF. If you use the empircal CDF array as "X" and the edges as "Y", then this is a PPF. In [12]: x = np.random.uniform(0, 10, size=1000) In [13]: hist, edges = np.histogram(x) In [15]: hist Out[15]: array([100, 101, 104, 108, 80, 111, 96, 88, 108, 104]) In [16]: edges Out[16]: array([ 3.53879571e-04, 1.00026423e+00, 2.00017458e+00, 3.00008493e+00, 3.99999528e+00, 4.99990563e+00, 5.99981598e+00, 6.99972633e+00, 7.99963668e+00, 8.99954704e+00, 9.99945739e+00]) In [18]: ecdf = np.hstack([0.0, hist.cumsum() / float(hist.sum())]) In [19]: ecdf Out[19]: array([ 0. , 0.1 , 0.201, 0.305, 0.413, 0.493, 0.604, 0.7 , 0.788, 0.896, 1. ]) In [20]: np.interp(np.linspace(0, 10), edges, ecdf) Out[20]: array([ 0. , 0.0203746 , 0.04078459, 0.06119459, 0.08160458, 0.10203472, 0.12264881, 0.14326291, 0.163877 , 0.18449109, 0.20522712, 0.22645351, 0.24767991, 0.2689063 , 0.29013269, 0.31160366, 0.33364646, 0.35568925, 0.37773204, 0.39977483, 0.41953158, 0.43585957, 0.45218756, 0.46851556, 0.48484355, 0.50433802, 0.52699311, 0.54964821, 0.5723033 , 0.59495839, 0.61577382, 0.63536742, 0.65496101, 0.6745546 , 0.6941482 , 0.71259664, 0.73055743, 0.74851823, 0.76647902, 0.78443982, 0.80567348, 0.82771627, 0.84975906, 0.87180185, 0.89384465, 0.91515087, 0.93637726, 0.95760365, 0.97883004, 1. ]) In [21]: np.interp(np.linspace(0, 1), ecdf, edges) Out[21]: array([ 3.53879571e-04, 2.04417216e-01, 4.08480553e-01, 6.12543890e-01, 8.16607227e-01, 1.02046852e+00, 1.22251143e+00, 1.42455434e+00, 1.62659724e+00, 1.82864015e+00, 2.02980301e+00, 2.22601775e+00, 2.42223250e+00, 2.61844725e+00, 2.81466200e+00, 3.01047705e+00, 3.19942458e+00, 3.38837211e+00, 3.57731965e+00, 3.76626718e+00, 3.95521472e+00, 4.19462069e+00, 4.44969986e+00, 4.70477903e+00, 4.95985820e+00, 5.15488346e+00, 5.33872431e+00, 5.52256515e+00, 5.70640600e+00, 5.89024684e+00, 6.08569264e+00, 6.29825861e+00, 6.51082459e+00, 6.72339057e+00, 6.93595654e+00, 7.16204944e+00, 7.39393960e+00, 7.62582975e+00, 7.85771991e+00, 8.07294833e+00, 8.26189586e+00, 8.45084340e+00, 8.63979093e+00, 8.82873846e+00, 9.01838365e+00, 9.21459840e+00, 9.41081315e+00, 9.60702789e+00, 9.80324264e+00, 9.99945739e+00]) -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco