code: def find_gamma_params(mode, std):
      
      if mode < 0:
          raise ValueError('input mode must be positive for gamma; got %e' % mode)
  
      var = std**2
      beta = (mode + np.sqrt(mode**2 + 4 * var)) / (2 * var)
      j = 2 * var / mode**2
      alpha = (j + 1 + np.sqrt(2 * j + 1)) / j
  
      if alpha <= 1:
          raise ValueError('couldn\'t compute self-consistent gamma parameters: '
                           'mode=%e std=%e alpha=%e beta=%e' % (mode, std, alpha, beta))
  
      return alpha, beta
docs: Given a modal value and a standard deviation, compute corresponding
      parameters for the gamma distribution.
  
      Intended to be used to replace normal distributions when the value must be
      positive and the uncertainty is comparable to the best value. Conversion
      equations determined from the relations given in the sample_gamma()
      docs.