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
def sample ( self ) : acceptance = 1 finish = 0 while ( acceptance < 0.234 or acceptance > 0.4 ) or finish == 0 : 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 ) 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 ) ) 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
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