mcmc_inversion.py
def do_mcmc(self, iterations, proposal_std):
        """ 
        Performs the MCMC sampling of the posterior probability
        """
        # Generate a model from the prior as an initial guess
        mcurr = self.generate_prior()
        # Remove the mean of the model
        mcurr -= np.mean(mcurr)

        # Array for storing the posterior distribution, posterior likelihoods,
        # and whether the trial was accepted/rejected
        posterior_distribution = np.zeros((self.mlen,iterations))
        posterior_likelihoods = np.zeros((iterations))
        accept_reject = np.zeros((100))
        for i in range(0, iterations):
            # Calculate the likelihood of the current model satisfying the data
            current_likelihood = self.evaluate_likelihood(mcurr)
            # Calculate the probability of the current model
            current_prob = current_likelihood

            # Calculate a new trial model
            mtrial = self.gaussian_mtrial(mcurr, proposal_std)
            # Remove the mean of the model
            mtrial -= np.mean(mtrial)
            # Check the prior probability of the proposed model
            prior_prob_trial = self.calculate_prior_uniform(mtrial)

            if not prior_prob_trial:
                # The trial is rejected
                accept_reject[i % 100] = False
                posterior_likelihoods[i] = current_likelihood
                posterior_distribution[:,i] = mcurr

                # Check for sensitivity test
                if not (i == 0) and (i % 100 == 0):
                    print(50*"-")
                    print("Sensitivity test for iteration number {}:".format(i))
                    proposal_std  = self.sensitivity_test(accept_reject, proposal_std)
                    print("The current likelihood is {}".format(
                        posterior_likelihoods[i]))
                    accept_reject = np.zeros((100))
                    print(50*"-")

                continue

            # Evaluate the likelihood of the trial model
            trial_likelihood = self.evaluate_likelihood(mtrial)
            # Calculate the probability of the trial model
            trial_prob = trial_likelihood
            
            # If the trial model has better or equal likelihood
            # than the current model, it becomes in the new
            # current model
            if trial_prob >= current_prob:
                #print("The likelihood is higher!")
                # The trial is accepted
                mcurr = mtrial
                accept_reject[i % 100] = True
                posterior_likelihoods[i] = trial_likelihood

            # Even if the likelihood is lower, we still may accept the
            # trial model. The probability that we do so is determined
            # the ratio of the two likelihoods
            else:
                # Likelihood ratio
                ratio = np.exp(trial_prob - current_prob)
                # Determine a random number between 0 and 1
                acceptance_thresh = np.random.rand(1)
                # If the likelihood ratio exceeds the acceptance threshold
                # We accept the new model
                if ratio > acceptance_thresh:
                    # The trial is accepted
                    mcurr = mtrial
                    accept_reject[i % 100] = True
                    posterior_likelihoods[i] = trial_likelihood
                else:
                    # The trial is rejected
                    #print("Proposal rejected due to likelihood")
                    accept_reject[i % 100] = False
                    posterior_likelihoods[i] = current_likelihood

            # Store the new model, whether the trial was accepted or rejected
            posterior_distribution[:,i] = mcurr

            if not (i == 0) and (i % 100 == 0):
                print(50*"-")
                print("Sensitivity test for iteration number {}:".format(i))
                proposal_std  = self.sensitivity_test(accept_reject, proposal_std)
                print("The current likelihood is {}".format(
                    posterior_likelihoods[i]))
                accept_reject = np.zeros((100))
                print(50*"-")

        return posterior_distribution, posterior_likelihoods
 
        Performs the MCMC sampling of the posterior probability
        
Performs the MCMC sampling of the posterior probability
Performs the MCMC sampling of the posterior probability
def do_mcmc(self, iterations, proposal_std):
        mcurr = self.generate_prior()
        mcurr -= np.mean(mcurr)
        posterior_distribution = np.zeros((self.mlen,iterations))
        posterior_likelihoods = np.zeros((iterations))
        accept_reject = np.zeros((100))
        for i in range(0, iterations):
            current_likelihood = self.evaluate_likelihood(mcurr)
            current_prob = current_likelihood
            mtrial = self.gaussian_mtrial(mcurr, proposal_std)
            mtrial -= np.mean(mtrial)
            prior_prob_trial = self.calculate_prior_uniform(mtrial)
            if not prior_prob_trial:
                accept_reject[i % 100] = False
                posterior_likelihoods[i] = current_likelihood
                posterior_distribution[:,i] = mcurr
                if not (i == 0) and (i % 100 == 0):
                    print(50*"-")
                    print("Sensitivity test for iteration number {}:".format(i))
                    proposal_std  = self.sensitivity_test(accept_reject, proposal_std)
                    print("The current likelihood is {}".format(
                        posterior_likelihoods[i]))
                    accept_reject = np.zeros((100))
                    print(50*"-")
                continue
            trial_likelihood = self.evaluate_likelihood(mtrial)
            trial_prob = trial_likelihood
            if trial_prob >= current_prob:
                mcurr = mtrial
                accept_reject[i % 100] = True
                posterior_likelihoods[i] = trial_likelihood
            else:
                ratio = np.exp(trial_prob - current_prob)
                acceptance_thresh = np.random.rand(1)
                if ratio > acceptance_thresh:
                    mcurr = mtrial
                    accept_reject[i % 100] = True
                    posterior_likelihoods[i] = trial_likelihood
                else:
                    accept_reject[i % 100] = False
                    posterior_likelihoods[i] = current_likelihood
            posterior_distribution[:,i] = mcurr
            if not (i == 0) and (i % 100 == 0):
                print(50*"-")
                print("Sensitivity test for iteration number {}:".format(i))
                proposal_std  = self.sensitivity_test(accept_reject, proposal_std)
                print("The current likelihood is {}".format(
                    posterior_likelihoods[i]))
                accept_reject = np.zeros((100))
                print(50*"-")
        return posterior_distribution, posterior_likelihoods
def do_mcmc ( self , iterations , proposal_std ) : mcurr = self . generate_prior ( ) mcurr -= np . mean ( mcurr ) posterior_distribution = np . zeros ( ( self . mlen , iterations ) ) posterior_likelihoods = np . zeros ( ( iterations ) ) accept_reject = np . zeros ( ( 100 ) ) for i in range ( 0 , iterations ) : current_likelihood = self . evaluate_likelihood ( mcurr ) current_prob = current_likelihood mtrial = self . gaussian_mtrial ( mcurr , proposal_std ) mtrial -= np . mean ( mtrial ) prior_prob_trial = self . calculate_prior_uniform ( mtrial ) if not prior_prob_trial : accept_reject [ i % 100 ] = False posterior_likelihoods [ i ] = current_likelihood posterior_distribution [ : , i ] = mcurr if not ( i == 0 ) and ( i % 100 == 0 ) : print ( 50 * "-" ) print ( "Sensitivity test for iteration number {}:" . format ( i ) ) proposal_std = self . sensitivity_test ( accept_reject , proposal_std ) print ( "The current likelihood is {}" . format ( posterior_likelihoods [ i ] ) ) accept_reject = np . zeros ( ( 100 ) ) print ( 50 * "-" ) continue trial_likelihood = self . evaluate_likelihood ( mtrial ) trial_prob = trial_likelihood if trial_prob >= current_prob : mcurr = mtrial accept_reject [ i % 100 ] = True posterior_likelihoods [ i ] = trial_likelihood else : ratio = np . exp ( trial_prob - current_prob ) acceptance_thresh = np . random . rand ( 1 ) if ratio > acceptance_thresh : mcurr = mtrial accept_reject [ i % 100 ] = True posterior_likelihoods [ i ] = trial_likelihood else : accept_reject [ i % 100 ] = False posterior_likelihoods [ i ] = current_likelihood posterior_distribution [ : , i ] = mcurr if not ( i == 0 ) and ( i % 100 == 0 ) : print ( 50 * "-" ) print ( "Sensitivity test for iteration number {}:" . format ( i ) ) proposal_std = self . sensitivity_test ( accept_reject , proposal_std ) print ( "The current likelihood is {}" . format ( posterior_likelihoods [ i ] ) ) accept_reject = np . zeros ( ( 100 ) ) print ( 50 * "-" ) return posterior_distribution , posterior_likelihoods
Performs the MCMC sampling of the posterior probability
Performs the MCMC sampling of the posterior probability
""" 
        Performs the MCMC sampling of the posterior probability
        """ # Generate a model from the prior as an initial guess # Remove the mean of the model # Array for storing the posterior distribution, posterior likelihoods, # and whether the trial was accepted/rejected # Calculate the likelihood of the current model satisfying the data # Calculate the probability of the current model # Calculate a new trial model # Remove the mean of the model # Check the prior probability of the proposed model # The trial is rejected # Check for sensitivity test # Evaluate the likelihood of the trial model # Calculate the probability of the trial model # If the trial model has better or equal likelihood # than the current model, it becomes in the new # current model #print("The likelihood is higher!") # The trial is accepted # Even if the likelihood is lower, we still may accept the # trial model. The probability that we do so is determined # the ratio of the two likelihoods # Likelihood ratio # Determine a random number between 0 and 1 # If the likelihood ratio exceeds the acceptance threshold # We accept the new model # The trial is accepted # The trial is rejected #print("Proposal rejected due to likelihood") # Store the new model, whether the trial was accepted or rejected
self
iterations
proposal_std