Performs random walk Metropolis-Hastings

        Parameters
        ----------
        scale : float
            Default starting scale

        nsims : int
            Number of simulations

        printer : Boolean
            Whether to print results or not

        method : str
            What type of MCMC

        cov_matrix: None or np.array
            Can optionally provide a covariance matrix for M-H.

def _mcmc_fit(self, scale=1.0, nsims=10000, printer=True, method="M-H", 
        cov_matrix=None, map_start=True, quiet_progress=False, **kwargs):
        """ Performs random walk Metropolis-Hastings

        Parameters
        ----------
        scale : float
            Default starting scale

        nsims : int
            Number of simulations

        printer : Boolean
            Whether to print results or not

        method : str
            What type of MCMC

        cov_matrix: None or np.array
            Can optionally provide a covariance matrix for M-H.
        """
        scale = 2.38/np.sqrt(self.z_no)

        # Get Mode and Inverse Hessian information
        if self.model_type in ['GPNARX', 'GPR', 'GP'] or map_start is True:
            y = self.fit(method='PML', printer=False)
            starting_values = y.z.get_z_values()

            # TODO: Bad use of a try/except - remove in future
            try:
                ses = np.abs(np.diag(y.ihessian))
                if len(ses[np.isnan(ses)]) != 0:
                    ses = np.ones(ses.shape[0])
                cov_matrix = np.zeros((len(ses), len(ses)))
                np.fill_diagonal(cov_matrix, ses)
            except:
                pass
        else:
            starting_values = self.latent_variables.get_z_starting_values()

        if method == "M-H":
            sampler = MetropolisHastings(self.neg_logposterior, scale, nsims, starting_values, 
                cov_matrix=cov_matrix, model_object=None, quiet_progress=quiet_progress)
            chain, mean_est, median_est, upper_95_est, lower_95_est = sampler.sample()
        else:
            raise Exception("Method not recognized!")

        if len(self.latent_variables.z_list) == 1:
            self.latent_variables.set_z_values(mean_est,'M-H',None,chain)
            mean_est = self.latent_variables.z_list[0].prior.transform(mean_est)
            median_est = self.latent_variables.z_list[0].prior.transform(median_est)
            upper_95_est = self.latent_variables.z_list[0].prior.transform(upper_95_est)
            lower_95_est = self.latent_variables.z_list[0].prior.transform(lower_95_est)        

        else:
            self.latent_variables.set_z_values(mean_est, 'M-H', None, chain)

            for k in range(len(chain)):
                mean_est[k] = self.latent_variables.z_list[k].prior.transform(mean_est[k])
                median_est[k] = self.latent_variables.z_list[k].prior.transform(median_est[k])
                upper_95_est[k] = self.latent_variables.z_list[k].prior.transform(upper_95_est[k])
                lower_95_est[k] = self.latent_variables.z_list[k].prior.transform(lower_95_est[k])        

        self.latent_variables.estimation_method = 'M-H'

        theta, Y, scores, states, states_var, X_names = self._categorize_model_output(mean_est)
    
        # Change this in future
        try:
            latent_variables_store = self.latent_variables.copy()
        except:
            latent_variables_store = self.latent_variables

        return MCMCResults(data_name=self.data_name,X_names=X_names,model_name=self.model_name,
            model_type=self.model_type, latent_variables=latent_variables_store,data=Y,index=self.index,
            multivariate_model=self.multivariate_model,objective_object=self.neg_logposterior, 
            method='Metropolis Hastings',samples=chain,mean_est=mean_est,median_est=median_est,lower_95_est=lower_95_est,
            upper_95_est=upper_95_est,signal=theta,scores=scores, z_hide=self._z_hide,max_lag=self.max_lag,
            states=states,states_var=states_var)