def plot_fit(self, intervals=True, **kwargs):
        """ Plots the fit of the Gaussian process model to the data

        Parameters
        ----------
        beta : np.array
            Contains untransformed starting values for latent variables
        
        intervals : Boolean
            Whether to plot uncertainty intervals or not
        
        Returns
        ----------
        None (plots the fit of the function)
        """
        import matplotlib.pyplot as plt
        import seaborn as sns

        figsize = kwargs.get('figsize',(10,7))

        date_index = self.index[self.max_lag:]
        expectation = self.expected_values(self.latent_variables.get_z_values())
        variance = self.variance_values(self.latent_variables.get_z_values())
        upper = expectation + 1.98*np.power(np.diag(variance),0.5)
        lower = expectation - 1.98*np.power(np.diag(variance),0.5)

        plt.figure(figsize=figsize) 

        plt.subplot(2, 2, 1)
        plt.title(self.data_name + " Raw")  
        plt.plot(date_index,self.data*self._norm_std + self._norm_mean,'k')

        plt.subplot(2, 2, 2)
        plt.title(self.data_name + " Raw and Expected") 
        plt.plot(date_index,self.data*self._norm_std + self._norm_mean,'k',alpha=0.2)
        plt.plot(date_index,self.expected_values(self.latent_variables.get_z_values())*self._norm_std + self._norm_mean,'b')

        plt.subplot(2, 2, 3)
        plt.title(self.data_name + " Raw and Expected (with intervals)")    

        if intervals == True:
            plt.fill_between(date_index, lower*self._norm_std + self._norm_mean, upper*self._norm_std + self._norm_mean, alpha=0.2)          
            
        plt.plot(date_index,self.data*self._norm_std + self._norm_mean,'k',alpha=0.2)
        plt.plot(date_index,self.expected_values(self.latent_variables.get_z_values())*self._norm_std + self._norm_mean,'b')

        plt.subplot(2, 2, 4)

        plt.title("Expected " + self.data_name + " (with intervals)")   

        if intervals == True:
            plt.fill_between(date_index, lower*self._norm_std + self._norm_mean, upper*self._norm_std + self._norm_mean, alpha=0.2)          
            
        plt.plot(date_index,self.expected_values(self.latent_variables.get_z_values())*self._norm_std + self._norm_mean,'b')

        plt.show()
Plots the fit of the Gaussian process model to the data

        Parameters
        ----------
        beta : np.array
            Contains untransformed starting values for latent variables
        
        intervals : Boolean
            Whether to plot uncertainty intervals or not
        
        Returns
        ----------
        None (plots the fit of the function)