code: def sample_gamma(alpha, beta, size):
      
  
      if alpha <= 0:
          raise ValueError('alpha must be positive; got %e' % alpha)
      if beta <= 0:
          raise ValueError('beta must be positive; got %e' % beta)
      return np.random.gamma(alpha, scale=1./beta, size=size)
docs: This is mostly about recording the conversion between Numpy/Scipy
      conventions and Wikipedia conventions. Some equations:
  
      mean = alpha / beta
      variance = alpha / beta**2
      mode = (alpha - 1) / beta [if alpha > 1; otherwise undefined]
      skewness = 2 / sqrt(alpha)