<jupyter_start><jupyter_text>Normal Distribution---The p.d.f. of a normal distribution is \begin{align*} & p(x|\mu,\sigma^2) = \frac1{\sqrt{2\pi\sigma^2}} \exp\left[-\frac{(x-\mu)^2}{2\sigma^2}\right], \\ & -\infty0, \\ &\mathrm{E}[X]=\mu,\quad \mathrm{Var}[X]=\sigma^2.\end{align*}<jupyter_code>import numpy as np
import scipy.stats as st
import scipy.optimize as opt
import pandas as pd
from IPython.display import display
from bokeh.io import show, output_notebook
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, HoverTool, Slider, Span
from bokeh.plotting import figure
output_notebook()
def normal_pdf_plot(doc):
    slider_mu = Slider(value=0.0, start=-7.0, end=7.0, step=0.1, title='\u03BC')
    slider_sigma = Slider(value=1.0, start=0.1, end=4.0, step=0.1, title='\u03C3')
    
    x = np.linspace(-7.0, 7.0, 1001)
    source = ColumnDataSource(data=dict(x=x, y=st.norm.pdf(x)))
    hover = HoverTool(tooltips=[('x', '@x{0.0000}'), ('density', '@y{0.0000}')])
    p = figure(plot_width=400, plot_height=300, x_range=(-7, 7), y_range=(0, 1.1),
               tools=[hover], toolbar_location=None)
    p.line('x', 'y', source=source, line_color='navy', line_width=2)
    p.xaxis.axis_label = 'Normal distribution'
    p.yaxis.axis_label = 'Probability density'
    p.xgrid.grid_line_color = p.ygrid.grid_line_color = p.outline_line_color = None   

    def update_pdf(attr, old, new):
        mu = slider_mu.value
        sigma = slider_sigma.value
        source.data['y'] = st.norm.pdf(x, loc=mu, scale=sigma)

    for params in [slider_mu, slider_sigma]:
        params.on_change('value', update_pdf)
    
    doc.add_root(column(row(slider_mu, slider_sigma, width=400), p))

# Your local jupyter server must be 'http://localhost:8888/'
show(normal_pdf_plot)<jupyter_output><empty_output><jupyter_text>Inverse Gamma Distribution and Student's $t$-distribution---The p.d.f. of an inverse gamma distribution is\begin{equation*} p(x|\alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{-(\alpha+1)}\exp\left(-\frac{\beta}{x}\right),\ x > 0,\ \alpha > 0,\ \beta > 0.\end{equation*}The p.d.f. of a (Student's) $t$-distribution is \begin{align*} & p(x|\nu,\mu,\sigma^2) = \frac{\Gamma\left(\frac{\nu+1}2\right)}{\Gamma\left(\frac{\nu}2\right)\sqrt{\pi\nu\sigma^2}} \left[1 + \frac{(x-\mu)^2}{\nu\sigma^2}\right]^{-\frac{\nu+1}2}, \\ & -\infty0,\ -\infty0, \\ & \mathrm{E}[X] = \mu,\ (\nu > 1),\ \mathrm{Var}[X] = \frac{\nu}{\nu-2}\sigma^2,\ (\nu>2).\end{align*}<jupyter_code>value_a = np.array([1.0, 3.0, 5.0, 5.0])
value_b = np.array([2.0, 2.0, 2.0, 1.0])
value_n = np.array([1.0, 2.0, 5.0])
styles = ['solid', 'dashed', 'dashdot', 'dotted']
colors = ['navy', 'firebrick', 'green', 'olive']
x1 = np.linspace(0, 2.3, 1001)
ig = figure(plot_width=400, plot_height=300, x_range=(0, 2.3), y_range=(0, 5),
            toolbar_location=None)
for index in range(value_a.size):
    a_i = value_a[index]
    b_i = value_b[index]
    plot_label = '\u03B1 = {0:3.1f}, \u03B2 = {1:3.1f}'.format(a_i, b_i)
    ig.line(x1, st.invgamma.pdf(x1, a_i, scale=b_i), line_color=colors[index],
            line_width=2, line_dash=styles[index], legend_label=plot_label)
ig.xaxis.axis_label = 'Inverse gamma distribution'
ig.yaxis.axis_label = 'Probability density'
ig.legend.click_policy = 'hide'
ig.legend.location = 'top_right'
ig.legend.border_line_color = ig.xgrid.grid_line_color = ig.ygrid.grid_line_color = ig.outline_line_color = None
x2 = np.linspace(-6.5, 6.5, 1001)
t = figure(plot_width=400, plot_height=300, x_range=(-6.5, 6.5), y_range=(0, 0.42),
           toolbar_location=None)
for index in range(value_n.size):
    n_i = value_n[index]
    plot_label = '\u03BD = {0:3.1f}'.format(n_i)
    t.line(x2, st.t.pdf(x2, n_i), line_color=colors[index], line_width=2,
           line_dash=styles[index], legend_label=plot_label)
t.line(x2, st.norm.pdf(x2), line_color=colors[-1], line_width=2, 
       line_dash=styles[-1], legend_label='\u03BD = \u221E')
t.xaxis.axis_label = 't distribution'
t.yaxis.axis_label = 'Probability density'
t.legend.click_policy = 'hide'
t.legend.location = 'top_right'
t.legend.border_line_color = t.xgrid.grid_line_color = t.ygrid.grid_line_color = t.outline_line_color = None
show(row(ig, t))<jupyter_output><empty_output><jupyter_text>Normal-Inverse-Gamma Prior---The natural conjugate prior distribution for $(\mu,\sigma^2)$ is\begin{equation*} \mu|\sigma^2 \sim \mathrm{Normal}\left(\mu_0,\frac{\sigma^2}{n_0}\right), \  \sigma^2 \sim \mathrm{Inv. Gamma}\left(\frac{\nu_0}2,\frac{\lambda_0}2\right).\end{equation*}The joint p.d.f. of the prior distribution is given by\begin{align*} p(\mu,\sigma^2) &= p(\mu|\sigma^2)p(\sigma^2), \\ p(\mu|\sigma^2) &= \sqrt{\frac{n_0}{2\sigma^2}}\exp\left[-\frac{n_0(\mu-\mu_0)^2}{2\sigma^2}\right] , \\ p(\sigma^2) &= \frac{\left(\frac{\lambda_0}2\right)^{\frac{\nu_0}2}}{\Gamma\left(\frac{\nu_0}2\right)}(\sigma^2)^{-\left(\frac{\nu_0}2+1\right)}\exp\left(-\frac{\lambda_0}{2\sigma^2}\right).\end{align*}The joint distribution of $(\mu, \sigma^2)$ is often called the normal-inverse-gamma distribution.  Likelihood---The likelihood of $(\mu, \sigma^2)$ is\begin{align*} p(D|\mu,\sigma^2) &= \prod_{i=1}^n p(x_i|\mu,\sigma^2) \\ &= \prod_{i=1}^n \frac1{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right] \\ &= (2\pi\sigma^2)^{-\frac{n}2}\exp\left[-\frac{\sum_{i=1}^n(x_i-\mu)^2}{2\sigma^2}\right].\end{align*}Since\begin{align*} \sum_{i=1}^n(x_i-\mu)^2 &= \sum_{i=1}^n(x_i-\bar x+\bar x-\mu)^2 \\ &= \sum_{i=1}^n\left\{(x_i-\bar x)^2 + 2(x_i-\bar x)(\bar x - \mu) + (\bar x - \mu)^2\right\} \\ &= \sum_{i=1}^n(x_i-\bar x)^2 + n(\bar x - \mu)^2,\end{align*}the likelihood is rewritten as\begin{equation*} p(D|\mu,\sigma^2) \propto (\sigma^2)^{-\frac{n}2}\exp\left[-\frac{\sum_{i=1}^n(x_i-\bar x)^2 + n(\bar x - \mu)^2}{2\sigma^2}\right].\end{equation*}  Joint Posterior Distribution---Applying Bayes' theorem, we have\begin{align*} & p(\mu,\sigma^2|D) \nonumber \\ &\propto p(D|\mu,\sigma^2)p(\mu|\sigma^2)p(\sigma^2) \\ &\propto (\sigma^2)^{-\frac{n}2}\exp\left[-\frac{\sum_{i=1}^n(x_i-\bar x)^2 + n(\bar x - \mu)^2}{2\sigma^2}\right] \nonumber \\ &\quad\times (\sigma^2)^{-\frac12}\exp\left[-\frac{n_0(\mu-\mu_0)^2}{2\sigma^2}\right] \times (\sigma^2)^{-\left(\frac{\nu_0}2+1\right)}\exp\left[-\frac{\lambda_0}{2\sigma^2}\right] \\ &\propto (\sigma^2)^{-\frac{n+\nu_0+3}2}\exp\Bigg[-\frac1{2\sigma^2}\bigg\{\sum_{i=1}^n(x_i-\bar x)^2 + n(\bar x - \mu)^2  \nonumber\\ &\qquad\qquad\qquad\qquad\qquad\qquad\qquad + n_0(\mu-\mu_0)^2 + \lambda_0\bigg\}\Bigg].\end{align*}By completing the square, we have\begin{align*} & n(\bar x - \mu)^2 + n_0(\mu-\mu_0)^2 \\ &= (n+n_0)\mu^2 - 2(n\bar x + n_0\mu_0)\mu + n\bar x^2 + n_0\mu_0^2 \\ &= (n+n_0)\left(\mu - \frac{n\bar x + n_0\mu_0}{n+n_0}\right)^2 + \frac{nn_0}{n+n_0}(\mu_0-\bar x)^2.\end{align*}Therefore the joint posterior distribution of $(\mu, \sigma^2)$ is derived as \begin{align*} p(\mu,\sigma^2|D) &\propto (\sigma^2)^{-\frac12}\exp\left[-\frac{n_*(\mu-\mu_*)^2}{2\sigma^2}\right] \\ &\quad\times (\sigma^2)^{-\left(\frac{\nu_*}2+1\right)}\exp\left(-\frac{\lambda_*}{2\sigma^2}\right),\end{align*}where\begin{align*} \mu_* &= \frac{n\bar x + n_0\mu_0}{n + n_0},\quad n_* = n + n_0,\quad \nu_* = n + \nu_0, \\ \lambda_* &= \sum_{i=1}^n(x_i-\bar x)^2 + \frac{nn_0}{n+n_0}(\mu_0-\bar x)^2 + \lambda_0.\end{align*}This is also a normal-inverse-gamma distribution.\begin{equation*} \mu|\sigma^2,D \sim \mathrm{Normal}\left(\mu_*,\frac{\sigma^2}{n_*}\right),\quad \sigma^2|D \sim \mathrm{Inv. Gamma}\left(\frac{\nu_*}2,\frac{\lambda_*}2\right).\end{equation*} Marginal Posterior Distributions---The marginal posterior distribution of $\mu$ is a Student's $t$-distribution:\begin{equation*} \mu|D\sim\mathrm{t}\left(\nu_*,\mu_*,\tau_*^2\right),\quad \tau_*^2 = \frac{\lambda_*}{\nu_*n_*}.\end{equation*}The marginal posterior distribution of $\sigma^2$ is an inverse gamma distribution:\begin{equation*} \sigma^2|D \sim \mathrm{Inv. Gamma}\left(\frac{\nu_*}2,\frac{\lambda_*}2\right).\end{equation*} HDPI of the inverse gamma distribution<jupyter_code>def invgamma_hpdi(ci0, alpha, beta, prob):
    def hpdi_conditions(v, a, b, p):
        eq1 = st.invgamma.cdf(v[1], a, scale=b) - st.invgamma.cdf(v[0], a, scale=b) - p
        eq2 = st.invgamma.pdf(v[1], a, scale=b) - st.invgamma.pdf(v[0], a, scale=b)
        return np.hstack((eq1, eq2))
    return opt.root(hpdi_conditions, ci0, args=(alpha, beta, prob)).x<jupyter_output><empty_output><jupyter_text>Posterior statistics of the parameters in the normal distribution<jupyter_code>def gaussian_stats(data, hyper_param, prob):
    mu0 = hyper_param['mu0']
    n0 = hyper_param['n0']
    nu0 = hyper_param['nu0']
    lam0 = hyper_param['lam0']
    n = data.size
    mean_data = data.mean()
    ssd_data = n * data.var()
    n_star = n + n0
    mu_star = (n * mean_data + n0 * mu0) / n_star
    nu_star = n + nu0
    lam_star = ssd_data + n * n0 / n_star * (mu0 - mean_data)**2 + lam0
    tau_star = np.sqrt(lam_star / nu_star / n_star)
    sd_mu = st.t.std(nu_star, loc=mu_star, scale=tau_star)
    ci_mu = st.t.interval(prob, nu_star, loc=mu_star, scale=tau_star)
    mean_sigma2 = st.invgamma.mean(0.5*nu_star, scale=0.5*lam_star)
    mode_sigma2 = lam_star / (nu_star + 2.0)
    median_sigma2 = st.invgamma.median(0.5*nu_star, scale=0.5*lam_star)
    sd_sigma2 = st.invgamma.std(0.5*nu_star, scale=0.5*lam_star)
    ci_sigma2 = st.invgamma.interval(prob, 0.5*nu_star, scale=0.5*lam_star)
    hpdi_sigma2 = invgamma_hpdi(ci_sigma2, 0.5*nu_star, 0.5*lam_star, prob)
    stats_mu = np.hstack((mu_star, mu_star, mu_star, sd_mu, ci_mu, ci_mu))
    stats_sigma2 = np.hstack((mean_sigma2, median_sigma2, mode_sigma2,
                              sd_sigma2, ci_sigma2, hpdi_sigma2))
    stats = np.vstack((stats_mu, stats_sigma2))
    stats_string = ['mean', 'median', 'mode', 'sd', 'ci (lower)', 'ci (upper)', 'hpdi (lower)', 'hpdi (upper)']
    param_string = ['$\\mu$', '$\\sigma^2$']
    results = pd.DataFrame(stats, index=param_string, columns=stats_string)
    post_param = dict(
        mu_star = mu_star,
        tau_star = tau_star,
        nu_star = nu_star,
        lam_star = lam_star
    )
    return results, post_param<jupyter_output><empty_output><jupyter_text>Plotting the posterior distribution of the parameters in the normal distribution<jupyter_code>def normal_posterior_plot(hyper_param, post_param, bounds):
    mu0 = hyper_param['mu0']
    n0 = hyper_param['n0']
    nu0 = hyper_param['nu0']
    lam0 = hyper_param['lam0']
    tau0 = np.sqrt(lam0 / nu0 / n0)
    mu_star = post_param['mu_star']
    tau_star = post_param['tau_star']
    nu_star = post_param['nu_star']
    lam_star = post_param['lam_star']

    # posterior distribution of mu
    x1 = np.linspace(bounds[0, 0], bounds[0, 1], 1001)
    source1 = ColumnDataSource(
        data = dict(
            x = x1,
            posterior = st.t.pdf(x1, nu_star, loc=mu_star, scale=tau_star),
            prior = st.t.pdf(x1, nu0, loc=mu0, scale=tau0)
        )
    )
    hover1 = HoverTool(
        tooltips = [
            ('\u03BC', '@x'),
            ('posterior', '@posterior'),
            ('prior', '@prior')
        ]
    )
    p1 = figure(plot_width=400, plot_height=300, tools=[hover1], toolbar_location=None)
    p1.line('x', 'posterior', source=source1, line_color='navy', line_width=2, legend_label='Posterior')
    p1.line('x', 'prior', source=source1, line_color='firebrick', line_width=2, line_dash='dotted', legend_label='Prior')
    p1.xaxis.axis_label = '\u03BC'
    p1.yaxis.axis_label = 'Probability density'
    p1.legend.click_policy = 'hide'
    p1.legend.location = 'top_right'
    p1.legend.border_line_color = p1.xgrid.grid_line_color = p1.ygrid.grid_line_color = p1.outline_line_color = None

    # posterior distribution of sigma^2
    x2 = np.linspace(bounds[1, 0], bounds[1, 1], 1001)
    source2 = ColumnDataSource(
        data = dict(
            x = x2,
            posterior = st.invgamma.pdf(x2, 0.5*nu_star, scale=0.5*lam_star), 
            prior = st.invgamma.pdf(x2, 0.5*nu0, scale=0.5*lam0)
        )
    )
    hover2 = HoverTool(
        tooltips = [
            ('\u03C3**2', '@x'),
            ('posterior', '@posterior'),
            ('prior', '@prior')
        ]
    )
    p2 = figure(plot_width=400, plot_height=300, tools=[hover2], toolbar_location=None)
    p2.line('x', 'posterior', source=source2, line_color='navy', line_width=2, legend_label='Posterior')
    p2.line('x', 'prior', source=source2, line_color='firebrick', line_width=2, line_dash='dotted', legend_label='Prior')
    p2.xaxis.axis_label = '\u03C3**2'
    p2.yaxis.axis_label = 'Probability density'
    p2.x_range.range_padding = 0
    p2.legend.click_policy = 'hide'
    p2.legend.location = 'top_right'
    p2.legend.border_line_color = p2.xgrid.grid_line_color = p2.ygrid.grid_line_color = p2.outline_line_color = None
    show(row(p1, p2))<jupyter_output><empty_output><jupyter_text>Application 1: Simulated Data---We use artificial data generated from the Normal distribution:$$ x_1,\dots,x_{50} \sim \mathrm{Normal}(1,2^2).$$The prior distribution is$$ \mu|\sigma^2 \sim \mathrm{Normal}\left(0,\frac{\sigma^2}{0.2}\right), \  \sigma^2 \sim \mathrm{Inv. Gamma}\left(\frac{5}2,\frac{7}2\right).$$The marginal prior distribution of $\mu$ is$$ \mu \sim \mathrm{t}(5,0,7).$$<jupyter_code>mu = 1.0
sigma = 2.0
n = 50
np.random.seed(99)
data = st.norm.rvs(loc=mu, scale=sigma, size=n)
hyper_param = dict(
    mu0 = 0.0,
    n0 = 0.2,
    nu0 = 5.0,
    lam0 = 7.0
)
prob = 0.95
results, post_param = gaussian_stats(data, hyper_param, prob)
display(results)
normal_posterior_plot(hyper_param, post_param, np.array([[-6.0, 6.0], [0.0, 10.0]]))<jupyter_output><empty_output><jupyter_text>Application 2: Weight of dimes--- Variables:+ mass - mass of dime in grams+ year - year the dime was minted Source:Dimes - Weights of a sample of dimes.https://vincentarelbundock.github.io/Rdatasets/datasets.html<jupyter_code>dime_data = pd.read_csv('Dimes.csv', index_col=0)
weight = dime_data['mass'].values
display(dime_data)
hyper_param_dime = dict(
    mu0 = 0.0,
    n0 = 0.00001,
    nu0 = 0.00001,
    lam0 = 0.00001
)
prob = 0.95
results_dime, post_param_dime = gaussian_stats(weight, hyper_param_dime, prob)
display(results_dime)
normal_posterior_plot(hyper_param_dime, post_param_dime, np.array([[2.22, 2.28], [0, 0.0015]]))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20145278450363197, 'fraction_numerical': 0.043375994465582844, 'mean_word_length': 3.9371584699453552, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 34, 'https://': 1, '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': '463204', 'license': 'MIT', 'n_tokens_mistral': 6012, 'n_tokens_neox': 5748, 'n_words': 972, 'path': 'Mathematics/Statistics/Statistics and Probability Python Notebooks/Concepts_Notebooks--probability_and_statistics-notebooks/notebook-a/example_normal.ipynb', 'repo_name': 'okara83/Becoming-a-Data-Scientist', 'size': '14455'}