docs: def marginalize(self, variables, inplace=True):
          """
          Modifies the distribution with marginalized values.
  
          Parameters
          ----------
          variables: iterator over any hashable object.
                  List of variables over which marginalization is to be done.
  
          inplace: boolean
                  If inplace=True it will modify the distribution itself,
                  else would return a new distribution.
  
          Returns
          -------
          GaussianDistribution or None :
                  if inplace=True (default) returns None
                  if inplace=False return a new GaussianDistribution instance
  
          Examples
          --------
          >>> import numpy as np
          >>> from pgmpy.factors.distributions import GaussianDistribution as GD
          >>> dis = GD(variables=['x1', 'x2', 'x3'],
          ...          mean=[1, -3, 4],
          ...          cov=[[4, 2, -2],
          ...               [2, 5, -5],
          ...               [-2, -5, 8]]))
          >>> dis.variables
          ['x1', 'x2', 'x3']
          >>> dis.mean
          array([[ 1],
                 [-3],
                 [ 4]])
          >>> dis.covariance
          array([[ 4,  2, -2],
                 [ 2,  5, -5],
                 [-2, -5,  8]])
  
          >>> dis.marginalize(['x3'])
          dis.variables
          ['x1', 'x2']
          >>> dis.mean
          array([[ 1.],
                 [-3.]]))
          >>> dis.covariance
          array([[4., 2.],
                 [2., 5.]])
          """
          if not isinstance(variables, list):
              raise TypeError("variables: Expected type list or array-like,"
                              "got type {var_type}".format(
                                  var_type=type(variables)))
  
          phi = self if inplace else self.copy()
  
          index_to_keep = [self.variables.index(var) for var in self.variables
                           if var not in variables]
  
          phi.variables = [phi.variables[index] for index in index_to_keep]
          phi.mean = phi.mean[index_to_keep]
          phi.covariance = phi.covariance[np.ix_(index_to_keep, index_to_keep)]
          phi._precision_matrix = None
  
          if not inplace:
              return phi