0

I need to generate a truncated gamma distribution pdf curve and histogram in Python 3.2 on win7. 

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as sps

shape, scale = 2., 2. # mean and dispersion
counter =0
s = []
upper_bound = 4
lower_bound  = 0.5
while (counter <= 1000):
    t = np.random.gamma(shape, scale, 1)
    if (lower_bound <= t <= upper_bound) :   
       s.append(t)
       counter += 1 

 count, bins, ignored = plt.hist(s, 50, normed=True)

 // this part take s very long time
 y = bins**(shape-1)*(np.exp(-bins/scale) /
                  (sps.gamma(shape)*scale**shape))

 plt.plot(bins, y, linewidth=2, color='r')
 plt.show()

I find that the following code takes a very long time and the pop-up figure window becomes non-responding.

"y = bins**(shape-1)*(np.exp(-bins/scale)/(sps.gamma(shape)*scale**shape))"

If I remove the lower and upper bounds for the gamma distribution, it runs very fast.

Any help would be appreciated.

user3356689
  • 55
  • 1
  • 9

1 Answers1

1

In the line with np.random.gamma you should not need the size=1 argument, you just want a float, so just do:

t = np.random.gamma(shape, scale)

It is taking a long time at the moment because each t you generate is an array with 1 element, and your s is a large nested array

What you've done in the while loop is already truncating your distribution! Although a much quicker way would be to get the range you want afterwards, i.e. replace your whole while loop with:

t = np.random.gamma(shape, scale, 1000) # 1000 entries in a gamma distribution
s = filter( lambda x: upper_bound>x>lower_bound, t ) # get entries within bounds
Higgs
  • 550
  • 5
  • 18
  • it works. Why ? And, i need to generate a truncated gamma distribution, could you please tell me how to do that ? Thanks ! – user3356689 Mar 14 '14 at 02:23
  • as I said, if you include the `size=1` argument, it will give you say `[0.5]` instead of `0.5`, and your `s` becomes `[[0.5],[0.4],...]` instead of `[0.5,0.4]` which is what you want. you can always try to print/plot s with/without the `size=1` argument to see the difference. the [documentation](http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.gamma.html) for more details: did you copy the code from the example? – Higgs Mar 14 '14 at 11:20
  • I made some changes to the code from the example. BTW, how to generate truncated gamma distribution ? Thanks ! – user3356689 Mar 14 '14 at 12:59
  • see my edit for truncation, make sure you understand the code =) finally, welcome to SO, you just need to accept my answer to say thanks – Higgs Mar 14 '14 at 13:03