<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Somme de variables aléatoires
Step1:  Somme de deux v.a. discrètes
Step2:  Paramètre $p$ différent (TODO)
TODO
Step3:  Somme de deux v.a. suivant une loi de poisson
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
Step4:  Somme de deux v.a. suivant une loi normale (convolution de lois normales)
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)$.
Step5:  Somme d'une loi de poisson (loi discrète) et d'une loi normale (loi continue)
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.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
Explanation: Somme de variables aléatoires
End of explanation
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);
Explanation: Somme de deux v.a. discrètes: cas général
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
Somme de deux v.a. discrètes indépendantes
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
Somme de deux v.a. suivant une loi de Bernoulli
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).
C.f. https://math.stackexchange.com/questions/1153576/addition-of-two-binomial-distribution.
Somme de deux v.a. suivant une loi binomiale
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:
* $X_1$ donne le nombre de "succès" sur $n_1$ expériences de Bernoulli $\mathcal{B}(p)$ e.g. le nombre de fois que la bille est allée à droite dans une planche de Galton de $n_1$ lignes.
* $X_2$ donne le nombre de "succès" sur $n_2$ expériences de Bernoulli $\mathcal{B}(p)$ e.g. le nombre de fois que la bille est allée à droite dans une planche de Galton de $n_2$ lignes.
* Donc la somme de $X_1$ et $X_2$ donne le nombre de "succès" sur $n_1 + n_2$ expériences de Bernoulli $\mathcal{B}(p)$ e.g. le nombre de fois que la bille est allée à droite dans une planche de Galton de $n_1 + n_2$ lignes.
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).
C.f. https://math.stackexchange.com/questions/1153576/addition-of-two-binomial-distribution.
Paramètre $p$ identique
$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.
End of explanation
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))
Explanation: Paramètre $p$ différent (TODO)
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$
End of explanation
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);
Explanation: Somme de deux v.a. suivant une loi de poisson
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.
End of explanation
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);
Explanation: Somme de deux v.a. suivant une loi normale (convolution de lois normales)
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)$.
End of explanation
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)
# 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)
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
dist = PoissonGaussian(k, mu, sigma)
y4 = np.array([dist.pdf(x) for x in x_normal])
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
# 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');
#np.diff(y4_cum) / (x_normal[1] - x_normal[0])
Explanation: Somme d'une loi de poisson (loi discrète) et d'une loi normale (loi continue)
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.
https://math.stackexchange.com/questions/455624/sum-of-poisson-and-gaussian-random-variable
https://fr.mathworks.com/matlabcentral/answers/287934-adding-two-different-distributions-example-gaussian-and-poisson-distribution?requestedDomain=www.mathworks.com
HS mais intéressant quand même:
* https://stats.stackexchange.com/questions/86402/variance-of-the-sum-of-a-poisson-distributed-random-number-of-normally-distribu
* https://en.wikipedia.org/wiki/Compound_Poisson_distribution
End of explanation
<END_TASK>