How to plot gamma distribution with alpha and beta parameters in python

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.

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)