docs: def to_joint_gaussian(self):
          """
          Return an equivalent Joint Gaussian Distribution.
  
          Examples
          --------
  
          >>> import numpy as np
          >>> from pgmpy.factors.continuous import CanonicalDistribution
          >>> phi = CanonicalDistribution(['x1', 'x2'], np.array([[3, -2], [-2, 4]]),
                                    np.array([[5], [-1]]), 1)
          >>> jgd = phi.to_joint_gaussian()
          >>> jgd.variables
          ['x1', 'x2']
          >>> jgd.covariance
          array([[ 0.5  ,  0.25 ],
                 [ 0.25 ,  0.375]])
          >>> jgd.mean
          array([[ 2.25 ],
                 [ 0.875]])
  
          """
          covariance = np.linalg.inv(self.K)
          mean = np.dot(covariance, self.h)
  
          return GaussianDistribution(self.variables, mean, covariance)