Initializes the appoximate distibution for the model

        Parameters
        ----------
        phi : np.ndarray
            Latent variables

        start_diffuse: boolean
            Whether to start from diffuse values or not
        
        gaussian_latents: LatentVariables object
            Latent variables for the Gaussian approximation

        Returns
        ----------
        BBVI fit object

def initialize_approx_dist(self, phi, start_diffuse, gaussian_latents):
        """ Initializes the appoximate distibution for the model

        Parameters
        ----------
        phi : np.ndarray
            Latent variables

        start_diffuse: boolean
            Whether to start from diffuse values or not
        
        gaussian_latents: LatentVariables object
            Latent variables for the Gaussian approximation

        Returns
        ----------
        BBVI fit object
        """     

        # Starting values for approximate distribution
        for i in range(len(self.latent_variables.z_list)):
            approx_dist = self.latent_variables.z_list[i].q
            if isinstance(approx_dist, fam.Normal):
                self.latent_variables.z_list[i].q.mu0 = phi[i]
                self.latent_variables.z_list[i].q.sigma0 = np.exp(-3.0)

        q_list = [k.q for k in self.latent_variables.z_list]

        # Get starting values for states
        T, Z, R, Q = self._ss_matrices(phi)
        H, mu = self.family.approximating_model(phi, T, Z, R, Q, gaussian_latents.get_z_values(transformed=True)[0], self.data)
        a, V = self.smoothed_state(self.data, phi, H, mu)
        V[0][0][0] = V[0][0][-1] 

        for item in range(self.data_length):
            if start_diffuse is False:
                q_list.append(fam.Normal(a[0][item], np.sqrt(np.abs(V[0][0][item]))))
            else:
                q_list.append(fam.Normal(self.family.itransform(np.mean(self.data)), np.sqrt(np.abs(V[0][0][item]))))

        return q_list