code: def MatrixGamma(symbol, alpha, beta, scale_matrix):
      """
      Creates a random variable with Matrix Gamma Distribution.
  
      The density of the said distribution can be found at [1].
  
      Parameters
      ==========
  
      alpha: Positive Real number
          Shape Parameter
      beta: Positive Real number
          Scale Parameter
      scale_matrix: Positive definite real square matrix
          Scale Matrix
  
      Returns
      =======
  
      RandomSymbol
  
      Examples
      ========
  
      >>> from sympy.stats import density, MatrixGamma
      >>> from sympy import MatrixSymbol, symbols
      >>> a, b = symbols('a b', positive=True)
      >>> M = MatrixGamma('M', a, b, [[2, 1], [1, 2]])
      >>> X = MatrixSymbol('X', 2, 2)
      >>> density(M)(X).doit()
      exp(Trace(Matrix([
      [-2/3,  1/3],
      [ 1/3, -2/3]])*X)/b)*Determinant(X)**(a - 3/2)/(3**a*sqrt(pi)*b**(2*a)*gamma(a)*gamma(a - 1/2))
      >>> density(M)([[1, 0], [0, 1]]).doit()
      exp(-4/(3*b))/(3**a*sqrt(pi)*b**(2*a)*gamma(a)*gamma(a - 1/2))
  
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Matrix_gamma_distribution
  
      """
      if isinstance(scale_matrix, list):
          scale_matrix = ImmutableMatrix(scale_matrix)
      return rv(
          symbol, MatrixGammaDistribution, (alpha, beta, scale_matrix)
      )
  
imports: from math import prod
  from sympy.core.basic import Basic
  from sympy.core.numbers import pi
  from sympy.core.singleton import S
  from sympy.functions.elementary.exponential import exp
  from sympy.functions.special.gamma_functions import multigamma
  from sympy.core.sympify import sympify, _sympify
  from sympy.matrices import ImmutableMatrix, Inverse, Trace, Determinant, MatrixSymbol, MatrixBase, Transpose, MatrixSet, matrix2numpy
  from sympy.stats.rv import _value_check, RandomMatrixSymbol, NamedArgsMixin, PSpace, _symbol_converter, MatrixDomain, Distribution
  from sympy.external import import_module
label_desc: Use image node_id 2 for calling a global function with example usage: MatrixGamma(symbol, alpha, beta, scale_matrix) and returns: rv