<jupyter_start><jupyter_text>Markov Chain Monte Carlo for Continuous Variable Problem DefinitionData is for personnel change from last year to this year for 10 companies. Model is defined as follows:$$y_i | \mu \overset{iid}{\sim} N(\mu, 1)$$$$\mu \sim t(0, 1, 1)$$where yi represents personnel change for company i, and the distribution of yi given $\mu$ is a Normal distribution with mean $\mu$ and variance 1. Prior distribution of $\mu$ is a t distribution with location 0, scale parameter 1, and degrees of freedom 1 (also known as Cauchy's distribution).Model is not conjugate, thus posterior is not a standard form that we can conveniently sample. To get posterior samples, we will need to setup a Markov chain, whose stationery distribution is the posterior distribution we want.The posterior distribution has this form (derived in Lesson 2), this is our g($\mu$):$$p(\mu | y_1, ..., y_n) \propto \frac{exp[n(\overline{y}\mu - \mu^2/2)]}{1 + \mu^2}$$Because of numerical stability concerns, we will work with log posteriors:$$log(g(\mu)) = n(\overline{y}\mu - \mu^2/2) - log(1 + \mu^2)$$ Implementation<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
%matplotlib inline
def lg(mu, n, ybar):
    mu2 = np.power(mu, 2)
    return n * (ybar * mu - (mu2 / 2)) - np.log(1 + mu2)


def mh(n, ybar, n_iter, mu_init, cand_sd):
    mu_out = np.zeros((n_iter,))
    accpt = 0
    mu_now = mu_init
    lg_now = lg(mu_init, n, ybar)
    for i in range(n_iter):
        # draw candidate for parameter mu from proposal distribution
        mu_cand = stats.norm.rvs(loc=mu_now, scale=cand_sd, size=1)[0]
        # since normal distribution is symmetric, q distributions cancel out
        lg_cand = lg(mu_cand, n, ybar)
        alpha = np.exp([lg_cand - lg_now])
        # accept candidate with probability alpha
        u = stats.uniform.rvs(loc=0, scale=1, size=1)
        if u < alpha:
            mu_now = mu_cand
            accpt += 1
            lg_now = lg_cand
            
        mu_out[i] = mu_now
        
    return mu_out, accpt/n_iter
y = np.array([1.2, 1.4, -0.5, 0.3, 0.9, 2.3, 1.0, 0.1, 1.3, 1.9])
ybar = np.mean(y)
n = len(y)<jupyter_output><empty_output><jupyter_text>Visualizing data and prior distributionsPlot below superimposes the data (histogram), dots (actual data points) and sample mean (vertical dashed line), with the prior distribution (green). We can make the following inferences from the plot.* Discrepancy between our prior $\mu$ and what the data says $\mu$ should be. Prior suggests that most of the probability mass of $\mu$ is near 0, while the data suggests that most of it is near 1.* As a result, we expect the posterior to have most of its mass somewhere between the prior mean and the data mean.<jupyter_code>plt.xlim((-1.0, 3.0))
# histogram of y values
plt.hist(y, density=True, bins=5, label="y")
# actual data points along x axis
plt.scatter(y, np.zeros((n)), marker='o', color='r')
# sample mean
plt.axvline(ybar, 0, 1, color='r', linestyle='--', label="sample mean")
# prior (t distribution with df=1)
xs = np.arange(-1.0, 3.0, 0.1)
plt.plot(xs, stats.t.pdf(xs, df=1), color='g', label="prior")
plt.legend(loc="best")
_ = plt.show()<jupyter_output><empty_output><jupyter_text>Sampling from the posteriorWe will now run our Metropolis-Hastings sampler `mh` to generate samples from the posterior. Acceptance Rate too low<jupyter_code>np.random.seed(43)
post_mu, accpt_rate = mh(n, ybar, n_iter=1000, mu_init=0, cand_sd=3)
def traceplot(post_mu, accpt_rate):
    """ poor man's traceplot, in R they do it with the "coda" package """
    plt.figure(figsize=(10, 5))
    plt.title("acceptance rate: {:.1f}%".format(accpt_rate * 100))
    plt.xlabel("iterations")
    plt.plot(np.arange(len(post_mu)), post_mu)
    _ = plt.show()    

traceplot(post_mu, accpt_rate)<jupyter_output><empty_output><jupyter_text>The traceplot gives some basic feedback about whether the Markov Chain has reached the stationery distribution.In this case, it appears that our proposal step size was a little too large, because we have an acceptance rate of 14% which is below the recommended range of 23-50%. This shows up in the traceplot as well -- mu gets stuck at a single value, especially in the earlier iterations. In other words, our Markov chain is not moving around as much as we would like it to. We would like to increase the accceptance rate. Acceptance Rate too highTo increase the acceptance rate, we want to rerun the sampler with a different standard deviation. Usually, if we **increase the variability** in the distribution that creates the candidates, it will **decrease** the acceptance rate. So we need to use a smaller standard deviation to increase the acceptance rate.<jupyter_code>post_mu, accpt_rate = mh(n, ybar, n_iter=1000, mu_init=0, cand_sd=0.05)
traceplot(post_mu, accpt_rate)<jupyter_output><empty_output><jupyter_text>Our acceptance rate is now very high at 95%. Notice that we now see the opposite behavior. The step sizes are small and the chain wanders. It is going to take the chain a long time to explore the posterior distribution and converge to the stationery distribution. Acceptance Rate within acceptable limits<jupyter_code>post_mu, accpt_rate = mh(n, ybar, n_iter=1000, mu_init=0, cand_sd=0.9)
traceplot(post_mu, accpt_rate)<jupyter_output><empty_output><jupyter_text>Our acceptance rate is now within the good limits. Traceplot indicates that the sampler is exploring the posterior distribution, even though it occasionally gets stuck in a particular state. Initial Value too far from true valueWhat happens if we change our initial value far away from the true value.<jupyter_code>post_mu, accpt_rate = mh(n, ybar, n_iter=1000, mu_init=30, cand_sd=0.9)
traceplot(post_mu, accpt_rate)<jupyter_output><empty_output><jupyter_text>Here we can see the effect of using a crazy starting value. It took about a 100 iterations to find the stationery distribution. Once it hits the stationery distribution, it bounces around the stationery distribution. So if we discard the first 100 values, we have the posterior distribution we wanted. Visualizing the posterior distribution<jupyter_code>post_mu_keep = post_mu[100:]
plt.xlim((-1.0, 3.0))
plt.hist(post_mu_keep, density=True, label="posterior")
# prior (copied from previous plot)
xs = np.arange(-1.0, 3.0, 0.1)
plt.plot(xs, stats.t.pdf(xs, df=1), color='g', label="prior")
# sample mean (copied from previous plot)
plt.axvline(ybar, 0, 1, color='r', linestyle='--', label="sample mean")
plt.legend(loc="best")
_ = plt.show()
np.mean(post_mu_keep)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10434258142340169, 'fraction_numerical': 0.02080820265379976, 'mean_word_length': 4.477291494632535, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 29, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '199829', 'license': 'Apache-2.0', 'n_tokens_mistral': 2210, 'n_tokens_neox': 2075, 'n_words': 812, 'path': 'techniques-and-models/w02-04b-mcmc-demo-continuous.ipynb', 'repo_name': 'sujitpal/bayesian-stats-examples', 'size': '6632'}