<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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 -------- ['x1', 'x2', 'x3'] array([[ 1. ], [-4.5], [ 8.5]]) 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)