Sample from M-H algorithm

        Returns
        ----------
        chain : np.array
            Chains for each parameter

        mean_est : np.array
            Mean values for each parameter

        median_est : np.array
            Median values for each parameter

        upper_95_est : np.array
            Upper 95% credibility interval for each parameter

        lower_95_est : np.array
            Lower 95% credibility interval for each parameter

def sample(self):
        """ Sample from M-H algorithm

        Returns
        ----------
        chain : np.array
            Chains for each parameter

        mean_est : np.array
            Mean values for each parameter

        median_est : np.array
            Median values for each parameter

        upper_95_est : np.array
            Upper 95% credibility interval for each parameter

        lower_95_est : np.array
            Lower 95% credibility interval for each parameter           
        """     

        acceptance = 1
        finish = 0

        while (acceptance < 0.234 or acceptance > 0.4) or finish == 0:

            # If acceptance is in range, proceed to sample, else continue tuning
            if not (acceptance < 0.234 or acceptance > 0.4):
                finish = 1
                if not self.quiet_progress:
                    print("")
                    print("Tuning complete! Now sampling.")
                sims_to_do = self.nsims
            else:
                sims_to_do = int(self.nsims/2) # For acceptance rate tuning

            # Holds data on acceptance rates and uniform random numbers
            a_rate = np.zeros([sims_to_do,1])
            crit = np.random.rand(sims_to_do,1)
            post = multivariate_normal(np.zeros(self.param_no), self.cov_matrix)
            rnums = post.rvs()*self.scale
            
            for k in range(1,sims_to_do): 
                rnums = np.vstack((rnums,post.rvs()*self.scale))

            self.phi, a_rate = metropolis_sampler(sims_to_do, self.phi, self.posterior, 
                a_rate, rnums, crit)

            acceptance = a_rate.sum()/a_rate.shape[0]
            self.scale = self.tune_scale(acceptance,self.scale)
            if not self.quiet_progress:
                print("Acceptance rate of Metropolis-Hastings is " + str(acceptance))

        # Remove warm-up and thin
        self.phi = self.phi[int(self.nsims/2):,:][::self.thinning,:]

        chain = np.array([self.phi[i][0] for i in range(0, self.phi.shape[0])])

        for m in range(1, self.param_no):
            chain = np.vstack((chain, [self.phi[i][m] for i in range(0,self.phi.shape[0])]))

        if self.param_no == 1:
            chain = np.array([chain])

        mean_est = np.array([np.mean(np.array([self.phi[i][j] for i in range(0,self.phi.shape[0])])) for j in range(self.param_no)])
        median_est = np.array([np.median(np.array([self.phi[i][j] for i in range(0,self.phi.shape[0])])) for j in range(self.param_no)])
        upper_95_est = np.array([np.percentile(np.array([self.phi[i][j] for i in range(0,self.phi.shape[0])]), 95) for j in range(self.param_no)])
        lower_95_est = np.array([np.percentile(np.array([self.phi[i][j] for i in range(0,self.phi.shape[0])]), 5) for j in range(self.param_no)])

        return chain, mean_est, median_est, upper_95_est, lower_95_est