In [ ]:
%matplotlib inline
In [ ]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
Soit $X$ et $Y$ deux v.a. discrètes de lois respectives $\{(x_i,p_i) ; i \in I\}$ et $\{(y_j,q_j) ; j \in J\}$. La v.a. $Z=X+Y$ est aussi une v.a. discrète dont la loi de probabilité est définie par l'ensemble des valeurs possibles $\{(x_i + y_j) ; i \in I, j \in J\}$ et par les probabilités associées:
$$ P(Z=z_k) = \sum\left\{ \frac{P(X=x_i,Y=y_j)}{x_i+y_j} = z_k \right\} $$ce qui suppose de connaitre la loi du couple $(X,Y)$.
TODO: clarifier la notation précédente
TODO: exemple
Soit $X$ et $Y$ deux v.a. discrètes indépendantes de lois respectives $\{(x_i,p_i) ; i \in I\}$ et $\{(y_j,q_j) ; j \in J\}$. La v.a. $Z=X+Y$ est aussi une v.a. discrète dont la loi de probabilité est définie par l'ensemble des valeurs possibles $\{(x_i + y_j) ; i \in I, j \in J\}$ et par les probabilités associées:
$$ \begin{eqnarray} P(Z=z_k) & = \sum_{i \in I} P(X=x_i)P(Y=z_k-x_i) \\ & = \sum_{j \in J} P(Y=y_j)P(X=z_k-y_j) \end{eqnarray} $$TODO: clarifier la notation précédente
TODO: exemple
La somme de deux v.a. suivant une loi de Bernoulli est juste un cas particulier de la section suivante (somme de deux v.a. suivant une loi binomiale).
Soit $X_1 \sim \mathcal{B}(1, p_1)$ et $X_2 \sim \mathcal{B}(1, p_2)$ deux v.a. indépendantes suivant une loi de Bernoulli.
Si $p_1 = p_2$ alors la somme de ces deux v.a. suit une loi binomiale de paramètre $n = 1 + 1 = 2$ et $p$ : $X_1 + X_2 \sim \mathcal{B}(2, p)$. Exemple intuitif: imaginer une planche de Galton de $1$ rangée à laquelle on ajouterais une seconde rangée...
Si $p_1 \neq p_2$ alors la somme de ces deux v.a. ne suit plus une loi binomiale mais une loi poisson binomiale (loi de probabilité discrète de la somme d'épreuves de Bernoulli indépendantes de paramètre $p$ différents).
Soit $X_1 \sim \mathcal{B}(n_1, p_1)$ et $X_2 \sim \mathcal{B}(n_2, p_2)$ deux v.a. indépendantes suivant une loi binomiale.
Si $p_1 = p_2$ alors la somme de ces deux v.a. suit une loi binomiale de paramètre $n = n_1 + n_2$ et $p$ : $X_1 + X_2 \sim \mathcal{B}(n_1 + n_2, p)$. Exemple intuitif: imaginer une planche de Galton de $n_1$ rangées à laquelle on ajouterais $n_2$ rangées:
Si $p_1 \neq p_2$ alors la somme de ces deux v.a. ne suit plus une loi binomiale mais une loi poisson binomiale (loi de probabilité discrète de la somme d'épreuves de Bernoulli indépendantes de paramètre $p$ différents).
$X_1 \sim \mathcal{B}(n_1, p_1)$
$X_2 \sim \mathcal{B}(n_2, p_2)$
$X_1 + X_2 \sim \mathcal{B}(n_1 + n_2, p)$ si $X_1$ et $X_2$ sont indépendantes.
In [ ]:
p = 0.5
n1 = 5
n2 = 8
# Empirical distribution
num_samples = 1000000
x1 = np.random.binomial(n=n1, p=p, size=num_samples)
x2 = np.random.binomial(n=n2, p=p, size=num_samples)
x3 = np.random.binomial(n=n1+n2, p=p, size=num_samples)
x1x2 = x1 + x2
# Probability mass function
xmin = 0
xmax = n1 + n2
x = np.arange(xmin, xmax, 1)
dist1 = scipy.stats.binom(n=n1, p=p)
dist2 = scipy.stats.binom(n=n2, p=p)
dist3 = scipy.stats.binom(n=n1+n2, p=p)
y1 = dist1.pmf(x)
y2 = dist2.pmf(x)
y3 = dist3.pmf(x)
# Plots
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
ax.hist(x1, bins=list(range(n1+2)), label=r"$X1 \sim \mathcal{B}(n_1, p)$", alpha=0.5, normed=True, color="blue")
ax.hist(x2, bins=list(range(n2+2)), label=r"$X2 \sim \mathcal{B}(n_2, p)$", alpha=0.5, normed=True, color="red")
ax.hist(x3, bins=list(range(n1+n2+2)), label=r"$X3 \sim \mathcal{B}(n_1 + n_2, p)$", alpha=0.5, normed=True, color="green")
ax.hist(x1x2, bins=list(range(n1+n2+2)), label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
ax.plot(x, y1, 'b.', label=r"PMF $\mathcal{B}(n_1, p)$")
ax.plot(x, y2, 'r.', label=r"PMF $\mathcal{B}(n_2, p)$")
ax.plot(x, y3, 'g.', label=r"PMF $\mathcal{B}(n_1+n_2, p)$")
ax.legend(prop={'size': 18}, loc='best', fancybox=True, framealpha=0.5);
TODO:
$X_1 \sim \mathcal{B}(n, p_1)$
$X_2 \sim \mathcal{B}(n, p_2)$
$X_1 + X_2 \sim ...$ si $X_1$ et $X_2$ sont indépendantes.
$E(X_1 + X_2) = n p_1 + n p_2$
$Var(X_1 + X_2) = n p_1 (1 − p_1) + n p_2 (1 − p_2)$
The characteristic function is $(p_1 e^{it} + 1 − p_1)^n (p_2 e^{it} + 1 − p_2)^n$
In [ ]:
p1 = 0.5
p2 = 0.8
n = 10
# Empirical distribution
num_samples = 1000000
x1 = np.random.binomial(n=n, p=p1, size=num_samples)
x2 = np.random.binomial(n=n, p=p2, size=num_samples)
x3 = np.random.binomial(n=2*n, p=p1*p2, size=num_samples)
x1x2 = x1 + x2
# Probability mass function
xmin = 0
xmax = 2 * n
x = np.arange(xmin, xmax, 1)
dist1 = scipy.stats.binom(n=n, p=p1)
dist2 = scipy.stats.binom(n=n, p=p2)
y1 = dist1.pmf(x)
y2 = dist2.pmf(x)
#y3 = (p1 * e^{it} + 1. − p1)**n * (p2 * e^{it} + 1. − p2)**n
# Plots
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
ax.hist(x1, bins=list(range(n+2)), label=r"$X1 \sim \mathcal{B}(n, p_1)$", alpha=0.5, normed=True, color="blue")
ax.hist(x2, bins=list(range(n+2)), label=r"$X2 \sim \mathcal{B}(n, p_2)$", alpha=0.5, normed=True, color="red")
ax.hist(x1x2, bins=list(range(2*n+2)), label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
ax.plot(x, y1, 'b.', label=r"PMF $\mathcal{B}(n, p_1)$")
ax.plot(x, y2, 'r.', label=r"PMF $\mathcal{B}(n, p_2)$")
ax.legend(prop={'size': 18}, loc='best', fancybox=True, framealpha=0.5)
print("E(X_1 + X_2) = ", x1x2.mean())
print("n p_1 + n p_2 = ", n * p1 + n * p2)
print()
print("Var(X_1 + X_2) = ", x1x2.var())
print("n p_1 (1 − p_1) + n p_2 (1 − p_2) =", n * p1 * (1-p1) + n * p2 * (1-p2))
Si $X_1 \sim \mathcal{P}(\lambda_1)$ et $X_2 \sim \mathcal{P}(\lambda_2)$ sont indépendantes, alors $X_1 + X_2 \sim \mathcal{P}(\lambda_1 + \lambda_2)$.
C.f. https://fr.wikipedia.org/wiki/Loi_de_Poisson#Stabilit.C3.A9_de_la_loi_de_Poisson_par_la_somme.
In [ ]:
lambda1 = 3
lambda2 = 8
# Empirical distribution
num_samples = 1000000
x1 = np.random.poisson(lam=lambda1, size=num_samples)
x2 = np.random.poisson(lam=lambda2, size=num_samples)
x3 = np.random.poisson(lam=lambda1+lambda2, size=num_samples)
x1x2 = x1 + x2
# Probability mass function
xmin = 0
xmax = 30 # TODO
x = np.arange(xmin, xmax, 1)
dist1 = scipy.stats.poisson(lambda1)
dist2 = scipy.stats.poisson(lambda2)
dist3 = scipy.stats.poisson(lambda1+lambda2)
y1 = dist1.pmf(x)
y2 = dist2.pmf(x)
y3 = dist3.pmf(x)
# Plots
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
ax.hist(x1, bins=list(range(xmax+2)), label=r"$X1 \sim \mathcal{P}(\lambda_1)$", alpha=0.5, normed=True, color="blue")
ax.hist(x2, bins=list(range(xmax+2)), label=r"$X2 \sim \mathcal{P}(\lambda_2)$", alpha=0.5, normed=True, color="red")
ax.hist(x3, bins=list(range(xmax+2)), label=r"$X3 \sim \mathcal{P}(\lambda_1 + \lambda_2)$", alpha=0.5, normed=True, color="green")
ax.hist(x1x2, bins=list(range(xmax+2)), label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
ax.plot(x, y1, 'b.', label=r"PMF $\mathcal{P}(\lambda_1)$")
ax.plot(x, y2, 'r.', label=r"PMF $\mathcal{P}(\lambda_2)$")
ax.plot(x, y3, 'g.', label=r"PMF $\mathcal{P}(\lambda_1 + \lambda_2)$")
ax.legend(prop={'size': 18}, loc='best', fancybox=True, framealpha=0.5);
La convolution (somme) de deux lois normales indépendantes forme une loi normale.
Soit $X_1 \sim \mathcal{N}(\mu_1, \sigma_1)$ et $X_2 \sim \mathcal{N}(\mu_2, \sigma_2)$ deux v.a. indépendantes suivant une loi Normale.
$X_1 + X_2 \sim \mathcal{N}\left(\mu_1 + \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2}\right)$.
In [ ]:
mu_1 = 1
mu_2 = 15
sigma_1 =2.5
sigma_2 = 1.5
# Empirical distribution
num_samples = 1000000
x1 = np.random.normal(loc=mu_1, scale=sigma_1, size=num_samples)
x2 = np.random.normal(loc=mu_2, scale=sigma_2, size=num_samples)
x3 = np.random.normal(loc=mu_1+mu_2, scale=math.sqrt(sigma_1**2 + sigma_2**2), size=num_samples)
x1x2 = x1 + x2
# Probability mass function
xmin = -10
xmax = 30
x = np.arange(xmin, xmax, 0.01)
dist1 = scipy.stats.norm(loc=mu_1, scale=sigma_1)
dist2 = scipy.stats.norm(loc=mu_2, scale=sigma_2)
dist3 = scipy.stats.norm(loc=mu_1+mu_2, scale=math.sqrt(sigma_1**2 + sigma_2**2))
y1 = dist1.pdf(x)
y2 = dist2.pdf(x)
y3 = dist3.pdf(x)
# Plots
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
bins = np.arange(xmin, xmax, 0.5)
ax.hist(x1, bins=bins, label=r"$X1 \sim \mathcal{N}(\mu_1, \sigma_1)$", alpha=0.5, normed=True, color="blue")
ax.hist(x2, bins=bins, label=r"$X2 \sim \mathcal{N}(\mu_2, \sigma_2)$", alpha=0.5, normed=True, color="red")
ax.hist(x3, bins=bins, label=r"$X3 \sim \mathcal{N}(\mu_1 + \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2})$", alpha=0.5, normed=True, color="green")
ax.hist(x1x2, bins=bins, label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
ax.plot(x, y1, 'b', label=r"PDF $\mathcal{N}(\mu_1, \sigma_1)$")
ax.plot(x, y2, 'r', label=r"PDF $\mathcal{N}(\mu_2, \sigma_2)$")
ax.plot(x, y3, 'g', label=r"PDF $\mathcal{N}(\mu_1 + \mu_2, \sqrt{\sigma_1^2 + \sigma_2^2})$")
ax.legend(prop={'size': 18}, loc='best', fancybox=True, framealpha=0.5);
Soit $X \sim \mathcal{P}(\lambda)$, $Y \sim \mathcal{N}(\mu, \sigma)$ deux v.a. indépendantes suivant respectivement une loi de poisson et une loi normale et $Z$ la somme de ces deux v.a.: $Z = X + Y$.
$X_1 + X_2 \sim \dots$ TODO
TODO $$PDF_Z(z) = \sum_{x=0}^{x=\infty} PMF_X(x) PDF_Y(z-x)$$
$E(Z) = $
$V(Z) = $
If X ~ Poisson(lambda), Y ~ N(mu,sigma^2), X, Y independent and Z=X+Y, then the cdf of Z is given by
$$P(Z \leq z) = \sum_{k=0}^{k=\infty} P(X=k) P(Y \leq z-k)$$avec
$$P(X=k) = \frac{\lambda^k}{k!} e^{-\lambda}$$et
$$P(Y \leq z-k) = 0.5 * \left(1 + \text{erf}\left(\frac{z - k - \mu}{\sqrt{(2 * \sigma^2)}}\right)\right)$$(erf: error function)
If needed, you can get the pdf of Z by differentiating the sum with respect to z.
In [ ]:
k = 2 # lambda
mu = 5
sigma = 0.5
# Empirical distribution
num_samples = 10000000
x1 = np.random.poisson(lam=k, size=num_samples)
x2 = np.random.normal(loc=mu, scale=sigma, size=num_samples)
x1x2 = x1 + x2
# Probability mass function
dist1 = scipy.stats.poisson(k)
dist2 = scipy.stats.norm(loc=mu, scale=sigma)
xmin = -30
xmax = 30
x_poisson = np.arange(0, xmax, 1)
x_normal = np.arange(xmin, xmax, 0.01)
y1 = dist1.pmf(x_poisson)
y2 = dist2.pdf(x_normal)
In [ ]:
# Plots
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
bins1 = np.arange(xmin, xmax, 1)
bins2 = np.arange(xmin, xmax, 0.1)
ax.hist(x1, bins=list(range(xmax+2)), label=r"$X1 \sim \mathcal{P}(\lambda)$", alpha=0.5, normed=True, color="blue")
ax.hist(x2, bins=bins2, label=r"$X2 \sim \mathcal{N}(\mu, \sigma)$", alpha=0.5, normed=True, color="red")
ax.hist(x1x2, bins=bins2, label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
ax.plot(x_poisson, y1, 'b.', label=r"PMF $\mathcal{P}(\lambda)$")
ax.plot(x_normal, y2, 'r', label=r"PDF $\mathcal{N}(\mu, \sigma)$")
ax.legend(prop={'size': 18}, loc='best', fancybox=True, framealpha=0.5)
In [ ]:
class PoissonGaussian:
def __init__(self, lambda_, mu, sigma):
self.lambda_ = lambda_
self.mu = mu
self.sigma = sigma
def pdf(self, x):
pdf = 0.
norm_dist = scipy.stats.norm(self.mu, self.sigma)
poisson_dist = scipy.stats.poisson(self.lambda_)
x_poisson = 0
while poisson_dist.cdf(x_poisson) < 0.999: # iterate over the X r.v. (Poisson distribution)
pdf += poisson_dist.pmf(x_poisson) * norm_dist.pdf(x-x_poisson)
x_poisson += 1
return pdf
In [ ]:
dist = PoissonGaussian(k, mu, sigma)
y4 = np.array([dist.pdf(x) for x in x_normal])
In [ ]:
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
ax.hist(x1x2, bins=bins2, label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
#ax.plot(x_normal, y4_cum, 'g'); # Print the CDF
ax.plot(x_normal, y4, 'g'); # Print the PDF
In [ ]:
# Approximation by a Normal distribution
print("X1+X2 mean =", x1x2.mean())
print("mu + lambda =", mu + k)
print()
print("X1+X2 std =", x1x2.std())
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(12, 12))
ax.hist(x1x2, bins=bins2, label=r"$X1 + X2$", alpha=0.5, normed=True, color="black", histtype="step", linewidth=2)
dist3 = scipy.stats.norm(loc=x1x2.mean(), scale=x1x2.std())
y3 = dist3.pdf(x_normal)
ax.plot(x_normal, y3, 'g');
In [ ]:
#np.diff(y4_cum) / (x_normal[1] - x_normal[0])