<jupyter_start><jupyter_text>Likelihood-Free Markov Chain Monte Carlo <NAME> In this tutorial We will write our own Markov Chain Monte Carlo algorithm with a Metropolis-Hastings sampler, then write a likelihood-free version. 1) MCMC with likelihoodsSometimes it's feasible to write down a [likelihood](https://en.wikipedia.org/wiki/Likelihood_function). In this example, we can use the $\chi^2$ to measure the goodness-of-fit, and use the [Metropolis-Hastings](https://en.wikipedia.org/wiki/Metropolis–Hastings_algorithm) algorithm to sample from the posterior distribution:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

np.random.seed(42)

# Set the properties of the line that we'll be fitting: 
true_std = 3
true_intercept = 0.5

x = np.linspace(0, 10, 500)
x -= x.mean()
y = true_std * np.random.randn(len(x)) + true_intercept
yerr = true_std * np.ones_like(x)

plt.errorbar(x, y, yerr, fmt='.')
plt.plot(x, true_intercept + np.zeros_like(x), color='k')
plt.xlabel('x')
plt.ylabel('y')
plt.show()<jupyter_output><empty_output><jupyter_text>We define a simple linear model to describe the data, which has parameters $\theta = \{m, b\}$.<jupyter_code>def linear_model(theta, x): 
    """
    Simple linear model. 
    """
    m, b = theta
    return m * x + b


def chi2(theta, x, y, yerr, model): 
    """
    Compute the \chi^2 by comparing `model` evaluated at `theta`
    to `y` at each value of `x`. 
    """
    return np.sum((model(theta, x) - y)**2 / yerr**2)<jupyter_output><empty_output><jupyter_text>Let the proposal step simply draw from a Gaussian distribution:<jupyter_code>def proposal(theta, scale=1):
    """
    Generate proposal step, by adding a draw from a 
    Gaussian distribution to the initial step position
    """
    return theta + scale * np.random.randn(len(theta))<jupyter_output><empty_output><jupyter_text>We will decide whether or not to accept new proposal steps using the Metropolis-Hastings algorithm, as defined by [Ford (2005)](https://arxiv.org/abs/astro-ph/0305441):<jupyter_code>def metropolis_hastings(init_theta, x, y, yerr, acceptance, scale=0.05):
    """
    Metropolis-Hastings algorithm, a la Ford (2005). 
    
    1) Generate a proposal step
    2) Compare the chi^2 of the proposal step to the current step
    3) Draw a random number `u` on [0, 1]
    4) Compute alpha = min([exp(-0.5 * (chi2(new) - chi2(old)), 1])
    5) If u <= alpha, accept step, otherwise keep step
    """
    # Generate a proposal step: 
    proposed_theta = proposal(init_theta, scale=scale)

    # Compare chi^2 of proposed step to current step:
    chi2_init_step = chi2(init_theta, x, y, yerr, linear_model)
    chi2_proposed_step = chi2(proposed_theta, x, y, yerr, linear_model)
    relative_likelihood = np.exp(-0.5 * (chi2_proposed_step - chi2_init_step))
    
    alpha = np.min([relative_likelihood, 1])

    # If U(0, 1) <= alpha, accept the step: 
    if np.random.rand() <= alpha: 
        return proposed_theta, acceptance + 1
    else: 
        return init_theta, acceptance
    

def sampler(x, y, yerr, init_theta, n_steps, scale=0.05):
    """
    Markov Chain Monte Carlo sampler. 
    """
    current_theta = np.copy(init_theta)
    # Allocate memory for samples: 
    samples = np.zeros((n_steps, len(init_theta)))
    samples[0, :] = init_theta
    acceptance = 0
    
    for i in range(1, n_steps):
        # Run the M-H algorithm to determine next step:
        current_theta, acceptance = metropolis_hastings(current_theta, 
                                                        x, y, yerr, acceptance, 
                                                        scale=scale)
        # Record the result: 
        samples[i, :] = current_theta
        
    # Compute the final acceptance rate
    acceptance_rate = acceptance / n_steps
    
    return samples, acceptance_rate<jupyter_output><empty_output><jupyter_text>Make an initial guess and let the MCMC algorithm find the correct solution:<jupyter_code>init_parameters = [0, 0.5]  # slope, intercept
n_steps = 50000

# This tweakable parameter determines how
# far new steps should be taken away from 
# previous steps. Increase `scale` to 
# decrease your acceptance rate: 
scale = 0.08

samples, acceptance_rate = sampler(x, y, yerr, init_parameters, 
                                   n_steps, scale=scale)<jupyter_output><empty_output><jupyter_text>What is the acceptance rate? Ideally this should be near 45%:<jupyter_code>print(acceptance_rate)<jupyter_output><empty_output><jupyter_text>Let's plot a few random draws from the posterior probability distribution functions for each parameter:<jupyter_code>for i in range(100): 
    random_step = np.random.randint(0, samples.shape[0])
    random_theta = samples[random_step, :]
    plt.plot(x, linear_model(random_theta, x), alpha=0.05, color='k')
plt.errorbar(x, y, yerr, fmt='.')
plt.xlabel('x')
plt.ylabel('y')<jupyter_output><empty_output><jupyter_text>You can see that the uncertainty in the measurements is being reflected by uncertainty in the slope and intercept parameters. We normally see the posterior probability distribution functions displayed in a "corner" plot like the one below:<jupyter_code>from corner import corner

burned_in_samples = samples[1000:]

corner(burned_in_samples, labels=['m', 'b'], truths=[0, true_intercept]);<jupyter_output><empty_output><jupyter_text>Note that the posterior distributions are consistent within the uncertainties for each parameter! We are accurately measuring each parameter and their uncertainties. In this example the likelihood was easy to compute, but there are cases when it's expensive or impossible to compute the likelihood. When that's the case, you can use... 2) *Likelihood-Free* MCMCFollowing the algorithm outlined by [Marjoram et al. 2003](http://www.pnas.org/content/100/26/15324).Likelihood-free MCMC bypasses the step within the M-H algorithm where you compare the likelihood of the observations given the parameters. Instead, you: 1. propose new parameters $\theta^\prime$2. generate a simulated dataset $y^\prime \equiv D^\prime$ from the new parameters3. measure the "distance" $\rho(D, D^\prime)$ between the simulated and observed datasets4. if the distance is less than some tolerance $h$, accept the step, return to (1)In this example, we'll be using a model that generates simulated datasets $D^\prime$ from a normal distribution, and computes the Euclidian *distance* between the means and standard deviations of the observations $D$ and the simulated dataset $D^\prime$, $$ \rho(D, D^\prime) = \sqrt{\left(\bar{y_\mathrm{obs}} - \bar{y^\prime}\right)^2 + \left(\mathrm{std}(y_\mathrm{obs}) - \mathrm{std}(y^\prime)\right)^2},$$and accept the step if $\rho(D, D^\prime) \leq h_j$ for $h \in \{1, 0.5, 0.2\}$.<jupyter_code>def D_prime(theta): 
    """
    Simulate a dataset by drawing from a random normal distribution
    """
    mean, std = theta
    return mean + std * np.random.randn(len(x))

def rho(theta, x, y, yerr):
    """
    Distance metric
    """
    # Generate a simulated dataset for comparison with observations
    simulated = D_prime(theta)
    
    # Compute distance between observations and simulation
    distance = np.sqrt((simulated.mean() - y.mean())**2 + (simulated.std() - y.std())**2)
    return distance

    
def metropolis_hastings_marjoram(init_theta, x, y, yerr, acceptance, scale, h):
    """
    Metropolis-Hastings algorithm, a la Marjoram et al. (2003)
    """
    # Generate a proposal step: 
    proposed_theta = proposal(init_theta, scale=scale)

    # Compute distance of proposed step
    rho_proposed_step = rho(proposed_theta, x, y, yerr)
    
    # If distance is less than h, accept step
    if rho_proposed_step <= h: 
        return proposed_theta, acceptance + 1

    else: 
        return init_theta, acceptance

def sampler_marjoram(x, y, yerr, init_theta, n_steps, scale, h):
    """
    Markov Chain Monte Carlo sampler without likelihoods 
    """
    current_theta = np.copy(init_theta)
    # Allocate memory for samples: 
    samples = np.zeros((n_steps, len(init_theta)))
    samples[0, :] = init_theta
    acceptance = 0
    
    for i in range(1, n_steps):
        # Run the M-H algorithm to determine next step:
        current_theta, acceptance = metropolis_hastings_marjoram(current_theta, 
                                                                 x, y, yerr, acceptance, 
                                                                 scale, h)
        # Record the result: 
        samples[i, :] = current_theta
        
    # Compute the final acceptance rate
    acceptance_rate = acceptance / n_steps
    
    return samples, acceptance_rate<jupyter_output><empty_output><jupyter_text>We'll loop over several values for $h$ and see how the posterior changes as $h \rightarrow 0$:<jupyter_code>init_parameters = [y.mean(), y.std()]  # mean, std
n_steps = 100000

# This tweakable parameter determines how
# far new steps should be taken away from 
# previous steps. Increase `scale` to 
# decrease your acceptance rate: 
scales = [1, 0.5, 0.05]

# This tweakable parameter sets the minimum
# distance required to accept a step in the MC
tolerances = [1, 0.5, 0.2]

burned_in_samples_h = []

for h, scale in zip(tolerances, scales):
    samples, acceptance_rate = sampler_marjoram(x, y, yerr, init_parameters, 
                                                n_steps, scale=scale, h=h)

    burned_in_samples_h.append(samples[n_steps//4:])

    print(f"h = {h:0.2f}, acceptance = {acceptance_rate:0.2f}")
bins = 35

histrange = [-0.5, 1.5]

fig, ax = plt.subplots(figsize=(4, 3))
ax.hist(burned_in_samples[:, 1], histtype='step', lw=2, 
         label='MCMC', density=True, range=histrange, bins=bins)

for s, h in zip(burned_in_samples_h, tolerances):
    ax.hist(s[:, 0], histtype='step', lw=2, bins=bins,
             label=f'LF-MCMC: h={h:.2f}', density=True, range=histrange)
ax.legend(loc=(1.01, 0.2))
for s in ['right', 'top']:
    ax.spines[s].set_visible(False)
    
ax.set(xlabel='Mean (intercept)', ylabel='PDF')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11012257809410834, 'fraction_numerical': 0.017002767892447607, 'mean_word_length': 3.484485815602837, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 42, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '360112', 'license': 'MIT', 'n_tokens_mistral': 3315, 'n_tokens_neox': 3109, 'n_words': 998, 'path': 'mcmc_likelihood_free_interact.ipynb', 'repo_name': 'bmorris3/mcmc_interact', 'size': '10122'}