21

I want to plot a gamma distribution with alpha = 29 (the scale) and beta = 3 (the size). In other words, I want to plot the pdf for Gamma(29,3). How do I do this if according to the documentation, the python gamma function only has parameters a and x and the size parameter doesn't exist?

I thought loc was beta, but I think it's actually offset, so the code below is wrong...

import numpy as np
import scipy.stats as stats 
from matplotlib import pyplot as plt

x = np.linspace (0, 100, 200) 
y1 = stats.gamma.pdf(x, a=29, loc=3) #a is alpha, loc is beta???
plt.plot(x, y1, "y-", label=(r'$\alpha=29, \beta=3$')) 


plt.ylim([0,0.08])
plt.xlim([0,150])
plt.show()
14wml
  • 4,048
  • 11
  • 49
  • 97

3 Answers3

26

According to the documentation, you want to use the scale parameter (theta), but since you are defining beta, which is the inverse of theta, then you pass scale with the value of 1/beta, which in your example would be 1/3 or 0.33333.

Therefore, try:

y1 = stats.gamma.pdf(x, a=29, scale=0.33333)
Scratch'N'Purr
  • 9,959
  • 2
  • 35
  • 51
  • 2
    so `a` is the shape in this case? – steven Dec 12 '18 at 17:35
  • 1
    @steven yes, it's the shape parameter – Scratch'N'Purr Dec 12 '18 at 19:06
  • 2
    It seems that scale is the same as beta, not the inverse. Take a look at these code snippets: `from math import e, gamma; fram scipy.stats import inv.gamma; a=.1; b=.5; x=4; (b**a)/gamma(a) * x**(-a-1) * e**(-b/x); invgamma.pdf(x, a, scale=b)` The latter two expressions equal the same! – Hielke Walinga Dec 18 '19 at 23:22
4

As @Hielke replied, as far as explained in scipy.stats 1.4.1 documentation it seems that the scalar parameter is equal to beta. Indeed, the function originally developped is :

gamma.pdf(x, a) = x^(a-1) * exp(-x) / gamma(a)

If one replaces x by a combination of the two optional parameters loc and scale as :

x = (y - loc) / scale

One should have :

gamma.pdf(x, a) = (y - loc)^(a-1) * exp( -(y - loc)/scale ) / (scale^(a-1) * gamma(a))

If you take loc = 0 then you recognized the expression of the Gamma distribution as usually defined. You multiply by the inverse of scale and you can conclude that scale = beta in this function and loc is an offset.

Actually I have tried to detail the documentation explaination :

Specifically, gamma.pdf(x, a, loc, scale) is identically equivalent to gamma.pdf(y, a) / scale with y = (x - loc) / scale.

eidal
  • 156
  • 2
  • 8
  • This seems to be misleading: beta is a rate, not scale parameter. As mentioned in the accepted answer, you want an inverse of the scale to get rate (beta). – Karaszka Jun 02 '20 at 11:17
  • 1
    The beta as used in the [Wiki article](https://en.wikipedia.org/wiki/Inverse-gamma_distribution) is equivalent to the scale parameter in `scipy.stats.invgamma`. Just do a quick plot of the mathematical expression for the pdf against values returned by `invgamma.pdf` to convince yourself. – marnix Feb 09 '21 at 17:51
1

This is not strictly answer for question, but this comes up when searching equivalent for R qgamma with alpha and beta parameters in Python. So just a note:

R:
qgamma(0.025, 5, 41.3)

Python:
from scipy.stats import gamma
gamma.ppf(0.025, 5, scale=1/41.3)
Napuu
  • 834
  • 8
  • 10