docs: def to_joint_gaussian(self):
          """
          The linear Gaussian Bayesian Networks are an alternative
          representation for the class of multivariate Gaussian distributions.
          This method returns an equivalent joint Gaussian distribution.
  
          Returns
          -------
          GaussianDistribution: An equivalent joint Gaussian
                                     distribution for the network.
  
          Reference
          ---------
          Section 7.2, Example 7.3,
          Probabilistic Graphical Models, Principles and Techniques
  
          Examples
          --------
          >>> from pgmpy.models import LinearGaussianBayesianNetwork
          >>> from pgmpy.factors.continuous import LinearGaussianCPD
          >>> model = LinearGaussianBayesianNetwork([('x1', 'x2'), ('x2', 'x3')])
          >>> cpd1 = LinearGaussianCPD('x1', [1], 4)
          >>> cpd2 = LinearGaussianCPD('x2', [-5, 0.5], 4, ['x1'])
          >>> cpd3 = LinearGaussianCPD('x3', [4, -1], 3, ['x2'])
          >>> model.add_cpds(cpd1, cpd2, cpd3)
          >>> jgd = model.to_joint_gaussian()
          >>> jgd.variables
          ['x1', 'x2', 'x3']
          >>> jgd.mean
          array([[ 1. ],
                 [-4.5],
                 [ 8.5]])
          >>> jgd.covariance
          array([[ 4.,  2., -2.],
                 [ 2.,  5., -5.],
                 [-2., -5.,  8.]])
  
          """
          variables = nx.topological_sort(self)
          mean = np.zeros(len(variables))
          covariance = np.zeros((len(variables), len(variables)))
  
          for node_idx in range(len(variables)):
              cpd = self.get_cpds(variables[node_idx])
              mean[node_idx] = sum([coeff * mean[variables.index(parent)] for
                                    coeff, parent in zip(cpd.beta_vector, cpd.evidence)]) + cpd.beta_0
              covariance[node_idx, node_idx] = sum(
                  [coeff * coeff * covariance[variables.index(parent), variables.index(parent)]
                   for coeff, parent in zip(cpd.beta_vector, cpd.evidence)]) + cpd.variance
  
          for node_i_idx in range(len(variables)):
              for node_j_idx in range(len(variables)):
                  if covariance[node_j_idx, node_i_idx] != 0:
                      covariance[node_i_idx, node_j_idx] = covariance[node_j_idx, node_i_idx]
                  else:
                      cpd_j = self.get_cpds(variables[node_j_idx])
                      covariance[node_i_idx, node_j_idx] = sum(
                          [coeff * covariance[node_i_idx, variables.index(parent)]
                           for coeff, parent in zip(cpd_j.beta_vector, cpd_j.evidence)])
  
          return GaussianDistribution(variables, mean, covariance)