I'm working through the Coursera course Bayesian Statistics: From Concept to Data Analysis by the University of California, Santa Cruz and there is a worksheet that requires me to use R or Excel to find the answer to this:

Let $Y ∼ Gamma(2, 1/3)$. Find $P(0.5 < Y < 1.5)$. The answer is given as 0.078.

I would like to calculate this using Python. I have tried

from scipy import stats
stats.gamma.cdf(1.5,1/3,scale=2) - stats.gamma.cdf(0.5,1/3,scale=2)

which returns 0.197. I've also tried switching the 2 and the 1/3.

Note I am using

Python 3.5.3 |Anaconda 2.5.0 (64-bit)| (default, May 15 2017, 10:43:23) [MSC v.1900 64 bit (AMD64)] on win32

What am I doing wrong?

  • $\begingroup$ Beware the important distinction between the gamma function and the gamma distribution.. Your title appears to conflate the two. The gamma cdf is the regularized incomplete gamma function; the plain "gamma function" appears as a normalizing constant in the gamma density. $\endgroup$
    – Glen_b
    Jun 8 '17 at 0:50
  • $\begingroup$ @Glen_b thanks, I guess I want Gamma distribution. Any idea how to get the desired result in Python? The docs are pretty unclear as to what the parameters are $\endgroup$
    – Dan
    Jun 8 '17 at 0:55
  • 1
    $\begingroup$ There are different ways to parameterize the gamma distribution. See en.wikipedia.org/wiki/Gamma_distribution. I'm sure the worksheet's usage does not match up with Python's default. $\endgroup$
    – Cliff AB
    Jun 8 '17 at 0:56
  • $\begingroup$ The question is using a different parameterization from the function you're calling. The following two posts should address all of your issues: stats.stackexchange.com/questions/85974/… and stats.stackexchange.com/questions/80833/… $\endgroup$
    – Glen_b
    Jun 8 '17 at 0:57
  • $\begingroup$ OK so after trying many combinations based on the wikipedia article Cliff posted, this is the parameterisation in Python that matches R: stats.gamma.cdf(1.5,2,scale=3) - stats.gamma.cdf(0.5,2,scale=3) $\endgroup$
    – Dan
    Jun 8 '17 at 1:08