# import sympy, Gamma, density, Symbol, pprint
from sympy.stats import Gamma, density
from sympy import Symbol, pprint
  
k = Symbol("k", positive = True)
theta = Symbol("theta", positive = True)
z = Symbol("z")
  
# using sympy.stats.Gamma() method
X = Gamma("x", k, theta)
gamVar = density((X)(z))
  
pprint(gamVar)

                 -z  
                -----
     -k  k - 1  theta
theta  *z     *e     
---------------------
       Gamma(k)  

# import sympy, Gamma, density, Symbol, pprint
from sympy.stats import Gamma, density
from sympy import Symbol, pprint
  
z = Symbol("z")
  
# using sympy.stats.Gamma() method
X = Gamma("x", 1 / 3, 45)
gamVar = density((X)(z))
  
pprint(gamVar)

           -z     
           ---    
   3 ____   45    
   \/ 75 *e       
------------------
    2/3           
15*z   *Gamma(1/3)