<jupyter_start><jupyter_text>Bayesian StatisticsBayes rule$$p(\theta| D) = \frac{p(D|\theta)p(\theta)}{p(D)}$$However, in the samples that we will soon see it is far more convenient to write the above as$$p(\theta| D) \propto p(D|\theta)p(\theta)$$since $p(D)$ is a constant **that does not depend on $\theta$**.However, we Bayesian analysis is far more involved in solving the following integral:$$p(D^*|D) = \int p(D^*|\theta)p(\theta|D) d\theta$$We denote $D^*$ to be test data, $D$ to be the train set and $\theta$ to be a hidden parameter. This can often be thought of as doing model averaging if $p(\theta|D)$ is thought of as a weighting function. Coin Tossing ExampleSuppose we have a biased coin where the true probability of landing heads is 0.6. We assume however, (for some odd reason that it is biased towards landing tails). A suitable prior for this would be a beta distribution which has the following functional form:$$p(\theta) = \frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha + \beta)}\theta^{\alpha - 1}(1-\theta)^{\beta - 1}$$and the distribution looks like (https://en.wikipedia.org/wiki/Beta_distribution):![beta](./images/Beta_distribution_pdf.svg.png)Since we believe it is left skewed we will choose $\alpha < \beta$. Let $\alpha = 2$, $\beta = 3$ and hence, **this forms our prior distribution**. Assume that after 10 coin tosses we observe 7 heads and 3 tails. Each (individual) coin toss is distributed as a bernoulli distribution, conditional on $\theta$, the probability of achieving a heads. At this point a frequentist method (I don't mean this in a derogative fashion) would state that the probability of receiving a heads is 0.7.The likelihood of the data can be states as follows:$$ p(y_1, \cdots, y_{10}| \theta) = \prod_{i=1}^{10} p(y_i|\theta)$$ since each draw is independent **given $\theta$** where, $p(y_i|\theta) = \theta^{1(y_i=1)}(1-\theta)^{1(y_i=0)}$. Hence in our example, the likelihood $ p(y_1, \cdots, y_{10}| \theta) = \theta^7(1-\theta)^3$. **Keep in mind that $\theta$ is considered fixed as far as the likelihood is concerned**, hence it is a bernoullii distribution, **not a beta distribution**. The posterior\begin{align}p(\theta | y) \propto & \quad p(y|\theta) p(\theta) \\\propto & \quad \theta^7(1-\theta)^3 \theta^2(1-\theta)^3 = \theta^{9}(1-\theta)^6\end{align}The normalising constant is irrelevant as we can see that the posterior $p(\theta | y)$ is a beta distribution simply by looking at its functional form. In fact the normalising constant (also known as the partition function) is simply,$$\int_0^1 \theta^{9}(1-\theta)^6 d\theta$$ since all probability distributions integrate out to one. In this case this turns out to be $\frac{\Gamma(9)\Gamma(6)}{\Gamma(15)}$. What is far more important to recognise than the normalising constant is that $\theta \sim Be(9, 6)$. The probability of $p(y^*=1|y)$\begin{align}p(y^*=1|y) &= \int p(y^*=1|\theta)p(\theta|y) d\theta\end{align}A frequentist method would have chosen the Maximum-A-Posteriori (MAP) estimate and plugged into $\theta$ without integrating it out. $$\begin{align}p(y^*=1|y) &= \int p(y^*=1|\theta)p(\theta|y) d\theta \\&= \int \theta \frac{\Gamma(15)}{\Gamma(9)\Gamma(6)} \theta^{9}(1-\theta)^6 d\theta\\&= \frac{\Gamma(15)}{\Gamma(9)\Gamma(6)} \int\theta^{10}(1-\theta)^6 d\theta\\&= \frac{\Gamma(15)}{\Gamma(9)\Gamma(6)} \frac{\Gamma(10)\Gamma(6)}{\Gamma(16)} \int\frac{\Gamma(16)}{\Gamma(10)\Gamma(6)}\theta^{10}(1-\theta)^6 d\theta \\&= \frac{\Gamma(15)}{\Gamma(9)\Gamma(6)} \frac{\Gamma(10)\Gamma(6)}{\Gamma(16)}  \quad 1 \\&= \frac{9}{15}\end{align}$$Note the rule $\Gamma(n) = (n-1)\Gamma(n-1)$ was used.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

import pymc3 as pm

from scipy.special import gamma

%matplotlib inline
alpha = 2
beta = 3
y_obs = np.zeros(10)
y_obs[:7] = 1
niter = 10000

with pm.Model() as model:
    θ = pm.Beta('θ', alpha=alpha, beta=beta)
    y = pm.Bernoulli('y', p=θ, observed=y_obs)
    
    trace = pm.sample(niter)
# norm_const = gamma(9)*gamma(6)/gamma(15)
theta = np.linspace(0,1,1000)
post_theta = (theta**9)*(1-theta)**6
C = np.sum(post_theta*(theta[1]-theta[0]))
post_theta = post_theta/C
plt.hist(trace['θ'][1000:],100, normed='true')
plt.plot(theta, post_theta)
plt.show()<jupyter_output><empty_output><jupyter_text>$$\int \theta p(\theta | y) d\theta \approx \sum_{\theta\sim p(\theta|y)} \theta$$<jupyter_code>np.mean(trace['θ'][1000:])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1770345596432553, 'fraction_numerical': 0.031438127090301006, 'mean_word_length': 4.841145833333333, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '95237', 'license': 'Apache-2.0', 'n_tokens_mistral': 1559, 'n_tokens_neox': 1469, 'n_words': 490, 'path': 'Bayesian_stats/Lesson 0 - Coin Toss.ipynb', 'repo_name': 'eagleanurag/deepschool.io', 'size': '4485'}