Commit c11cfe3f authored by Rémy Huet's avatar Rémy Huet 💻
Browse files

Questions 6 & 7

parent 748fdd99
%% Cell type:markdown id:de03950d tags:
# AOS1 - TP01
## Introduction to Bayesian Inference
### Imports
%% Cell type:code id:505a1f86 tags:
```
import numpy as np
import scipy as sp
import scipy.stats as spst
import matplotlib.pyplot as plt
import itertools
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
```
%% Cell type:markdown id:cee8b996 tags:
### Maximum Likehood estimation
#### Random sample generation
%% Cell type:code id:b757b943 tags:
```
n = 10
p = 0.7
reps = 200
sample = np.random.binomial(n, p, reps)
# Count results for binomial experience
b_unique, b_count = np.unique(sample, return_counts=True)
# Get results using the law
x = np.arange(11)
width = 0.35
fig, ax = plt.subplots()
rects1 = ax.bar(b_unique - width/2, b_count / reps, width, label='Sample')
rects2 = ax.bar(x + width/2, spst.binom.pmf(x, n, p), width, label='Theorical')
# Add some text for labels, title and custom x-axis tick labels, etc.
ax.set_ylabel('Proportion / probability')
ax.set_title('Sample and theorical repartition for binomial law')
ax.legend()
plt.show()
```
%% Cell type:code id:4ee353be tags:
```
l = 10
reps = 200
sample = np.random.poisson(l, reps)
# Count results
p_unique, p_count = np.unique(sample, return_counts=True)
# Get results using the law
# Range using quartiles
max_x = max(max(p_unique) + 1, spst.poisson.ppf(0.99, l))
min_x = min(min(p_unique) + 1, spst.poisson.ppf(0.01, l))
x = np.arange(min_x,
max_x)
width = 0.35
fig, ax = plt.subplots()
rects1 = ax.bar(p_unique - width/2, p_count / reps, width, label='Sample')
rects2 = ax.bar(x + width/2, spst.poisson.pmf(x, l), width, label='Theorical')
ax.set_ylabel('Proportion / probability')
ax.set_title('Sample and theorical repartition for poisson law')
ax.legend()
plt.show()
```
%% Cell type:code id:cd7ad873-5ce9-4d46-9c6f-1f2d539d3675 tags:
```
alpha = 2
beta = 4.2
reps = 600
sampling = 75
sample = np.random.beta(alpha, beta, reps)
max_x = max(max(sample), spst.beta.ppf(0.99, alpha, beta))
min_x = min(min(sample), spst.beta.ppf(0.01, alpha, beta))
x = np.linspace(min_x, max_x, sampling)
fig, ax = plt.subplots()
ax.hist(sample, bins=75, density=True, label='Sample')
ax.plot(x, spst.beta.pdf(x, alpha, beta), label = 'Theorical')
ax.set_ylabel('Density')
ax.set_title('Sample and theorical repartition for beta law')
ax.legend()
plt.show()
```
%% Cell type:code id:1d6d5ceb-5b2b-4749-b14d-c22b9ddaf76a tags:
```
k = 9
theta = 0.5
reps = 600
sampling = 75
sample = np.random.gamma(k, theta, reps)
max_x = max(max(sample), spst.gamma.ppf(0.99, k, scale=theta))
min_x = min(min(sample), spst.gamma.ppf(0.01, k, scale=theta))
x = np.linspace(min_x, max_x, sampling)
fig, ax = plt.subplots()
ax.hist(sample, bins=75, density=True, label='Sample')
ax.plot(x, spst.gamma.pdf(x, k, scale=theta), label = 'Theorical')
ax.set_ylabel('Count')
ax.set_title('Sample and theorical repartition for gamma law')
ax.legend()
plt.show()
```
%% Cell type:code id:28f78cb1-39c9-43ef-b771-f5acfddbeabd tags:
```
l = 1
reps = 600
sampling = 75
sample = np.random.exponential(l, reps)
max_x = max(max(sample), spst.expon.ppf(0.99, scale = 1/l))
min_x = min(min(sample), spst.expon.ppf(0.01, scale = 1/l))
x = np.linspace(min_x, max_x, sampling)
fig, ax = plt.subplots()
ax.hist(sample, bins=75, density=True, label='Sample')
ax.plot(x, spst.expon.pdf(x, scale = 1/l), label = 'Theorical')
ax.set_ylabel('Count')
ax.set_title('Sample and theorical repartition for exponential law')
ax.legend()
plt.show()
```
%% Cell type:code id:06cab6c5-b950-407c-bbb8-c420be891adb tags:
```
mu = -0.5
sigma2 = 2.3
reps = 600
sampling = 75
sample = np.random.normal(mu, sigma2, reps)
max_x = max(max(sample), spst.norm.ppf(0.99, mu, scale = sigma2))
min_x = min(min(sample), spst.norm.ppf(0.01, mu, scale = sigma2))
x = np.linspace(min_x, max_x, sampling)
fig, ax = plt.subplots()
ax.hist(sample, bins=75, density=True, label='Sample')
ax.plot(x, spst.norm.pdf(x, mu, scale = sigma2), label = 'Theorical')
ax.set_ylabel('Count')
ax.set_title('Sample and theorical repartition for normal law')
ax.legend()
plt.show()
```
%% Cell type:markdown id:25a8d33c-2450-41c8-8390-d2987f6e3070 tags:
#### Likelihood plot
We suppose that the tirages are independant
$$L(\theta | \underline{x}) = f(x_1, x_2, \dots, x_n; \theta) = \Pi_{i = 1}^{n}p(x_i | \theta)$$
%% Cell type:code id:178ca65c-229d-476a-b244-e56beec792d3 tags:
```
def loglike(law, params, samples):
like = 1
for sample in samples:
like = like * law.pdf(sample, params)
return np.log(like)
```
%% Cell type:markdown id:b9a5b7c6-1b2f-4485-a6cd-e5c498166cfe tags:
## Bayesian updating using conjugate priors
### Beta-binomial distribution
1) La densité de probabilité de la loi bêta est pour $x \in [0; 1], \alpha \gt 0, \beta \gt 0$ :
$$\frac{x^{\alpha - 1}(1 - x)^{\beta - 1}}{B(\alpha, \beta)}$$
où $B(x, y) = \int_0^1 t^{x - 1}(1 - t)^{y - 1}dt$
Son espérence est $\frac{\alpha}{\alpha + \beta}$, son mode $\frac{\alpha - 1}{\alpha + \beta - 2}$ si $\alpha \gt 1$ et $\beta \gt 1$ et sa variance $\frac{\alpha\beta}{(\alpha + \beta)^2(\alpha + \beta + 1)}$
%% Cell type:code id:2169d739-dbd4-4030-87f0-542bdbf49152 tags:
```
@interact(a = widgets.FloatSlider(min = 0.1, max = 10, step = 0.1), b = widgets.FloatSlider(min = 0.1, max = 10, step = 0.1))
def trace_beta(a, b):
x = np.linspace(0, 1, 100)
fig, ax = plt.subplots()
ax.plot(x, spst.beta.pdf(x, a, b))
ax.set_title('Repartition of beta function')
plt.show()
```
%% Cell type:code id:d2ee0d03 tags:
%% Cell type:markdown id:78974b3b tags:
Qu 7
X ~ $B(n, \theta)$ n connu et $\theta$ ~ beta($\alpha, \beta$)
On a $\theta \in [0;1]$
On peut calculer :
* La distribution jointe $\phi(x, t) = f(x | t)\pi_\theta(t)$
* La distribution marginale de X $f_X(x) = \int_{\Theta}\phi(x,t)dt$
* La distribution postérieure de $\theta$ sachant $X$ $\pi_\theta(t|x) = \frac{\phi(x,t)}{f_X(x)}$
%% Cell type:markdown id:918e429c tags:
On a :
$$\pi_\theta(t,\alpha,\beta) = \frac{t^{\alpha - 1}(1 - t)^{\beta - 1}}{B(\alpha,\beta)}$$
$$f_X(x|\theta) = \begin{pmatrix} n \\ x \end{pmatrix} \theta^x (1 - \theta)^{n - x} $$
La distribution jointe est :
$$\phi(x, t) = \begin{pmatrix} n \\ x \end{pmatrix} \frac{t^{\alpha + x - 1}(1 - t)^{\beta + n - x - 1}}{B(\alpha, \beta)}$$
La distribution marginale de X est :
$$f_X(x) = \int_0^1\phi(x, t)dt = \int_0^1\begin{pmatrix} n \\ x \end{pmatrix}\frac{t^{\alpha + x - 1}(1 - t)^{\beta + n - x - 1}}{B(\alpha, \beta)}dt = \frac{\begin{pmatrix} n \\ x \end{pmatrix}}{B(\alpha, \beta)} \int_0^1 t^{\alpha + x - 1}(1 - t)^{\beta + n - x - 1}$$
Soit $$f_X(x) = \frac{\begin{pmatrix} n \\ x \end{pmatrix}}{B(\alpha, \beta)} B(\alpha + x, \beta + n - x)$$
On en déduit la distribution postérieure
$$\pi_\theta(t | x) = \frac{t^{\alpha + x - 1}(1 - t)^{\beta + n - x - 1}}{B(\alpha + x, \beta + n - x)}$$
On retrouve la loi bêta corrigée par l'échantillon.
%% Cell type:code id:0af57523 tags:
```
def trace(a, b, n, x):
sample = np.linspace(0, 1, 100)
fig, ax = plt.subplots()
ax.plot(sample, spst.beta.pdf(sample, a, b), label='Prior distribution')
ax.plot(sample, spst.binom.pmf(x, n, sample), label="Likehood of theta")
ax.plot(sample, spst.beta.pdf(sample, a + x, b + n - x), label="Posterior distribution")
ax.legend()
ax.set_title("Prior distribution, likehood of theta and posterior distribution")
plt.show()
interact(trace,
a = widgets.FloatSlider(min = 0.1, max = 10, step = 0.1),
b = widgets.FloatSlider(min = 0.1, max = 10, step = 0.1),
n = widgets.IntSlider(min = 1, max = 20),
x = widgets.IntSlider(min = 0, max = 20)
)
```
%% Cell type:code id:b864cb67 tags:
```
```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment