```{r include=FALSE, echo=FALSE, eval=TRUE}
library(knitr)

options(width = 300)
knitr::opts_chunk$set(
  fig.width = 7, fig.height = 5, 
  out.width = "80%",
  fig.align = "center", 
  fig.path = "figures/04_discrete-distrib_",
  size = "tiny", 
  echo = TRUE, eval = TRUE, 
  warning = FALSE, message = FALSE, 
  results = TRUE, comment = "")
# knitr::asis_output("\\footnotesize")

dir.main <- "~/stat1"
dir.slides <- file.path(dir.main, "slides")
setwd(dir.slides)
```


## Distribution de probabilité discrète

On parle de ***distribution discrète*** pour désigner la distribution de probabilité de variables ne pouvant prendre que des valeurs discrètes (par opposition aux distributions continues). 

**Notes: ** 

- En probabilités la variable observée ($x$) représente généralement le nombre de succès d'une série d'observations. Elle prend donc généralement des valeurs naturelles ($x \in \mathbb{N}$). 
- La probabilité $\operatorname{P}(x)$ prend des valeurs réelles entre 0 et 1, mais sa distribution est discrète puisqu'elle n'est définie que pour des valeurs discrètes de $x$. On la représente généralement par une fonction en escalier. 

## Distribution géométrique

**Application:** temps d'attente jusqu'à la première réalisation d'un évènement au cours d'un schéma de Bernoulli. 

**Exemples:**

- Comptage du nombre de jets d'un dé ($x$) qui précèdent la première occurrence d'un 6 (l'occurrence n'est pas incluse dans le compte).
- Longueur d'une séquence d'ADN avant la première occurrence d'une cytosine 


## Fonction de masse de probabilité géométrique

La ***fonction de masse de probabilité*** (***Probability Mass Function***, ***PMF***) indique la probabilité d'observer un résultat élémentaire particulier. 

Pour la distribution géométrique, elle indique la probabilité d'observer exactement $x$ échecs avant le premier succès, au cours d'une série d'essais indépendants à probabilité de succès $p$. 

$$\operatorname{P}(X = x) = (1-p)^x \cdot p$$

Justification: 

- Probabilité d'échec pour un essai = $q = 1-p$ (évènements complémentaires)
- Schéma de Bernoulli $\rightarrow$ les essais sont indépendants $\rightarrow$ probabilité de la série est le produit des probabilités des résultats successifs. 
- On calcule donc le produit des probabilité des $x$ échecs initiaux et du succès au $(x+1)^\text{ème}$ essai.

*Note: la PMF est apparentée au concept de **densité** que nous verrons lorsque nous traiterons les distributions continues.*


## Fonction de masse de probabilité géométrique


```{r geometric_PMF, fig.width=9, fig.height=4, echo=FALSE, fig.cap="**Fonction de masse de la loi géométrique**. Gauche: ordonnée en échelle logarithmique. "}
################################################################
##
## Illustrations of the geometric distribution
##
## Running this script requires to first run the script config.R
##    source('http://pedagogix-tagc.univ-mrs.fr/courses/stat
##
## Single-command execution:
##   source(file.path(dir.R.files, 'geometric.R'))

# x11(width = 12, height = 15)

## Parameters
p <- 0.25      # Probability of success
n <- 30        # Waiting time (time of the first success)
x <- 0:n       # X values to be plotted

test.x <- 3    #  A particular X value where an arraw will be plotted


## Color specification (I cannot use colors for the book, but well for slides)
color.status <- "full"
if (color.status == "full") {
  plot.colors <- c('density'='blue',
              'CDF'='#00BB00',
              'dCDF'='purple'
              )
} else {
  plot.colors <- c('density'='black',
                   'CDF'='black',
                   'dCDF'='black'
  )
}

# ## Drawing parameters
par.ori <- par(no.readonly = TRUE)
par(mfrow=c(1,2))
par(las=1)
grid.color <- '#CCCCCC'

## grid sizes
x.grid <- 0:6*5
y.grid <- 0:20*0.05

## ##############################################################
## Waiting time before the first success in a Bernouli schema

################################################################
## Plot the probability mass function
plot(x - 0.5,
     dgeom(x,p),
     ylim=c(0,1),
     main="Geometric\nProbability mass function (PMF)",
     xlab='Waiting time',
     ylab='Probability: P(X=x)',
     panel.first=c(abline(v=x.grid,col=grid.color),
       abline(h=y.grid,col=grid.color)),
     col=plot.colors['density'],
     lwd=2, las=1,
     type='s'
     )
legend("topright", bty = "o", bg="white", 
       legend=paste(sep="", 'P(X=',test.x,') =', format(p)))

## Draw arrows to illustrate the relationship between density curve and probability
test.p <- dgeom(test.x,p)
arrows(test.x, 0,test.x,test.p,col='#888888',length=0.1,angle=25,code=2,lwd=2)
arrows(test.x, test.p,0,test.p,col='#888888',length=0.1,angle=25,code=2,lwd=2)
text(test.x,0, test.x,pos=4, cex=0.7)
text(0,test.p, round(test.p, digits=3),pos=1, cex=0.7)

################################################################
## Plot the probability mass function with logarithmic Y scale. Show an
## extended domain of X, to illustrate the capability of log scales to
## display small probabilities.
ext.x <- 0:(n+10)
plot(ext.x - 0.5,
     dgeom(ext.x,p),
     main=paste('PMF (log Y scale)'),
     xlab='Waiting time',
     ylab='Probability: P(X=x); log scale',
     ylim=c(min(dgeom(ext.x,p)),1),
     panel.first=grid(equilog=F,col=grid.color,lty='solid'),
     col=plot.colors['density'],
     lwd=2, las=1,
     type='s',
     log='y'
     )

## Draw arrows to illustrate the relationship between CDF and left tail probability
log.test.x = c(test.x, 29)
log.test.p <- round(dgeom(log.test.x, p),digits=c(3,5))
arrows(log.test.x, 1e-5,log.test.x,log.test.p,col='#888888',length=0.1,angle=25,code=2,lwd=2)
arrows(log.test.x, log.test.p,0,log.test.p,col='#888888',length=0.1,angle=25,code=2,lwd=2)
text(log.test.x,1e-5, log.test.x,pos=1, cex=0.7)
text(0,log.test.p, log.test.p,pos=1, cex=0.7)

par(mfrow=c(1,1))
```

## Queues de distribution et fonction de répartition

Les queues de la distribution sont les aires comprises sous la courbe de densité jusqu'à une certaine valeur (***queue gauche***) ou à partir d’une certaine valeur (***queue droite***). 

- La ***queue droite*** indique la probabilité d'obtenir un résultat ($X$) **inférieur ou égal** à une certaine valeur ($x$): $\operatorname{P}(X \le x$). 

    - **Définition: ** la ***fonction de répartition*** (***Cumulative Density Function***, ***CDF***) $\operatorname{P}(X \le x)$ indique la probabilité qu'une variable aléatoire $X$ prenne une valeur inférieure ou égale à une valeur donnée ($x$). Elle correspond à la queue gauche (en incluant la valeur $x$ considérée). 

- La **queue gauche** d'une distribution indique la probabilité d'observer un résultat **supérieur ou égal** à une certaine valeur: $\operatorname{P}(X \ge x$). 

    - Note: nous verrons ultérieurement l'utilisation de la queue droite de différentes distributions en tant que probabilité critique (***P value***), dans le cadre de tests d'enrichissement fonctionnel, sur-représentation de motifs, de détection de pics, $\ldots$  

## Queues de distribution et fonction de répartition

```{r goemetric_CDF, fig.width=7, fig.height=5.5, out.width="90%", echo=FALSE, fig.cap="**Queues et fonction de répartition de probabilité géométrique**. "}

par(mfrow = c(2,2))

################################################################
## Plot the probability mass function and highlight the left tail
left.tail.range <- 0:test.x
plot(x - 0.5,
     dgeom(x,p),
     ylim = c(0,1),
     main = paste('Left tail, X<=', test.x),
     xlab = 'Waiting time',
     ylab = 'Probability: P(X=x)',
     panel.first = c(abline(v = x.grid,
                            col = grid.color),
                     abline(h = y.grid,col = grid.color),
       rect(left.tail.range - 0.5, 0, left.tail.range + 0.5, dgeom(left.tail.range,p), density = NULL, angle = 45,
            col = '#BBBBBB', border = '#BBBBBB')
       ),
     col = plot.colors['density'],
     lwd = 2, las = 1,
     type = 's'
     )
legend('topright',
       legend = paste(
         'P(X<=', test.x, ') = ', 
         round(pgeom(test.x, p), digits = 3), sep = ''), 
       bty = 'o', bg  = '#FFFFFF')
axis(1, at = test.x, labels = test.x, 
     las = 1, col.ticks = '#888888')


################################################################
## Plot the probability mass function and highlight the right tail
right.tail.range <- test.x:n
plot(x - 0.5,
     dgeom(x,p),
     ylim = c(0,1),
     main = paste('Right tail, X>=', test.x),
     xlab = 'Waiting time',
     ylab = 'Probability: P(X=x)',
     panel.first = c(abline(v = x.grid, col = grid.color),
       abline(h = y.grid, col = grid.color),
       rect(right.tail.range - 0.5, 0, right.tail.range + 0.5,
            dgeom(right.tail.range,p), density = NULL, angle = 45,
            col = '#BBBBBB', border = '#BBBBBB')
       ),
     col = plot.colors['density'],
     lwd = 2, las = 1,
     type = 's'
     )
legend('topright',
       legend = paste('P(X>=',test.x,') = ', 
                      round(pgeom(test.x - 1, p, lower.tail = F), 
                            digits = 3), sep = ''), 
       bty = 'o', 
       bg = '#FFFFFF')
axis(1, at = test.x, labels = test.x, las = 1, 
     col.ticks = '#888888')

################################################################
## Plot the CDF function
plot(x - 0.5,
     pgeom(x,p),
     ylim = c(0,1),
     main = paste('Cumulative distribution function (CDF)'),
     xlab = 'Waiting time',
     ylab = 'CDF = P(X<=x)',
     panel.first = c(abline(v = x.grid, col = grid.color),
       abline(h = y.grid,col = grid.color)),
     col = plot.colors['CDF'],
     lwd = 2, las = 1,
     type = 's'
     )
legend('bottomright', 
       legend = paste('P(X<=',test.x,') = ', 
                    round(pgeom(test.x, p), digits = 3),  
                    sep = ''), bty = 'o', bg = '#FFFFFF')

## Draw arrows to illustrate the relationship between CDF and left tail probability
left.tail.p <- pgeom(test.x, p)
arrows(test.x, 0, test.x, left.tail.p, 
       col = '#888888', length = 0.1, 
       angle = 25, code = 2, lwd = 2)
arrows(test.x, left.tail.p, 0, left.tail.p, 
       col = '#888888', length = 0.1, angle = 25, code = 2, lwd = 2)
text(x = test.x, y = 0, labels = test.x, pos = 4, cex = 0.7)
text(0, left.tail.p, round(left.tail.p, digits = 3), pos = 3, cex = 0.7)


################################################################
## Plot the dCDF function
plot(x - 0.5,
     pgeom(x - 1, p, lower.tail = F), ## We use x-1 to obtain the inclusive dCDF rather than the default exclusive one
     ylim = c(0,1),
     main = paste('Decreasing CDF (dCDF)'),
     xlab = 'Waiting time',
     ylab = 'dCDF = P(X>=x)',
     panel.first = c(abline(v = x.grid, col = grid.color),
       abline(h = y.grid, col = grid.color)),
     col = plot.colors['dCDF'],
     lwd = 2, las = 1,
     type = 's')
legend('topright', 
       legend = paste(
         'P(X>=', test.x,') = ', 
         round(pgeom(test.x - 1, p, lower.tail = F), digits = 3), 
         sep = ''), bty = 'o', 
       bg = '#FFFFFF')

## Draw arrows to illustrate the relationship between dCDF and right tail probability
right.tail.p <- pgeom(test.x - 1, p, lower.tail = F)
arrows(test.x, 0, test.x, right.tail.p, col = '#888888', 
       length = 0.1, angle = 25, code = 2, lwd = 2)
arrows(test.x, right.tail.p, 0, right.tail.p, col = '#888888',
       length = 0.1, angle = 25, code = 2, lwd = 2)
text(test.x, 0, test.x,pos = 4, cex = 0.7)
text(0, right.tail.p, round(right.tail.p, digits = 3), 
     pos = 3, cex = 0.7)

## Show that the sum of CDF and dCDF differ from 1
## print(c(dgeom(test.x, p), left.tail.p, right.tail.p, right.tail.p + left.tail.p),digits=3)
par(par.ori)
par(mfrow = c(1,1))

```

## Distribution binomiale

La **distribution binomiale** indique la probabilité d'observer un certain nombre ($x$) de succès au cours d'une série de $n$ essais indépendants avec une probabilité de succès $p$ constante (schéma de Bernoulli).

**Fonction de masse de probabilité binomiale**

$$\operatorname{P}(X=x) = \binom{n}{x} \cdot  p^x  \cdot (1-p)^{n-x} = C_n^x  p^x  (1-p)^{n-x} = \frac{n!}{x!(n-x)!} p^x (1-p)^{n-x}$$

**Fonction de répartition binomiale**

$$\operatorname{P}(X \ge x)  =   \sum_{i=x}^{n}{P(X=i)} =  \sum_{i=x}^{n}{C_n^i p^i (1-p)^{n-i}}$$

**Propriétés**

- **Espérance** (nombre de succès attendus au hasard): $<X> = n \cdot p$
- **variance**: $\sigma^2 = n \cdot p \cdot (1 - p)$. 

    - Note: la variance de la binomiale est inférieure à sa moyenne. 
    
- **Ecart-type**: $\sigma = \sqrt{n \cdot p \cdot (1 - p)}$



## Distribution binomiale en $i$

La distribution binomiale peut prendre différentes formes selon les valeurs des paramètres (probabilité de succès $p$, et nombre d'essais $n$). 

Quand l'espérance ($p \cdot n$) est inférieure à 1, la distribution est monotone décroissante, et on la qualifie de ***distribution en forme de $i$***. 

```{r binomial_i-shaped, fig.width=9, fig.height=4, echo=FALSE, fig.cap="Distribution binomiale en forme de i.  "}


## Drawing parameters
par.ori <- par(no.readonly = TRUE)
par(mfrow=c(1,2))
grid.color <- '#CCCCCC'


## color specification (actually I don't use colors for the book
if (color.status == "full") {
  plot.col <- c('probability'='blue',
                'CDF'='darkgreen',
                'dCDF'='purple',
                'grid'='#CCCCCC'
                )
} else {
  plot.col <- c('probability'='black',
                'CDF'='black',
                'dCDF'='black',
                'grid'='#CCCCCC'
                )
}
line.types <- c(probability="solid",
                CDF="solid",
                dCDF="solid",
                grid="solid")

## Drawing parameters
par.ori <- par(no.readonly = TRUE)
par(las=1)

#### initialize parameters
p.values <- c(0.02,0.1,0.5,0.98) ## Success probabilities
n <- 20 ## Number of trials
x.display <- 20 ## Max value to display on the X axis
x <- c(0:n, n) ## Add an "n" at the end in order to have the last horizontal bar in the plot
x.shifted <- (0:(n+1))-0.5
x.grid <- x
y.grid <- seq(from=0,to=1,by=0.2)
  

#### Binomial functions
p <- 0.02
binomial.mean <- p*n

xmax <- 5 + 4*binomial.mean

plot(x.shifted,
     dbinom(x,n,p),
     main=paste("Binomial probability (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X=x)",
     col=plot.col["probability"],
     lty=line.types["probability"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)

plot(x.shifted,
     pbinom(x,n,p),
     main=paste("CDF (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X <= x)",
     col=plot.col["CDF"],
     lty=line.types["CDF"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)
par(par.ori)
par(mfrow=c(1,1))
```

## Distribution binomiale en cloche asymétrique

Quand la probabilité de succès relativement élevée mais inférieure à 0.5, la distribution prend une forme en cloche asymétrique. 

```{r binomial_distrib_bell-shaped, fig.width=9, fig.height=4, echo=FALSE, fig.cap="Distribution binomiale en forme de cloche asymétrique.  "}

## Drawing parameters
par(mfrow=c(1,2))
par(las=1)

################################################################
#
# Binomial series
#


#### Binomial functions
p <- 0.25
binomial.mean <- p*n

xmax <- 5 + 4*binomial.mean

plot(x.shifted,
     dbinom(x,n,p),
     main=paste("Binomial probability (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X=x)",
     col=plot.col["probability"],
     lty=line.types["probability"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)

plot(x.shifted,
     pbinom(x,n,p),
     main=paste("CDF (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X <= x)",
     col=plot.col["CDF"],
     lty=line.types["CDF"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)
par(par.ori)
par(mfrow=c(1,1))
```

## Distribution binomiale en cloche symétrique

Quand la probabilité de succès vaut 0.5, la distribution prend une forme en cloche symétrique. 

```{r binomial_distrib_symmetric-bell, fig.width=9, fig.height=4, echo=FALSE, fig.cap="Distribution binomiale en forme de cloche symétrique (p=0.5). "}

## Drawing parameters
par(mfrow=c(1,2))
par(las=1)

p <- 0.5
binomial.mean <- p*n

xmax <- 5 + 4*binomial.mean

plot(x.shifted,
     dbinom(x,n,p),
     main=paste("Binomial probability (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X=x)",
     col=plot.col["probability"],
     lty=line.types["probability"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)

plot(x.shifted,
     pbinom(x,n,p),
     main=paste("CDF (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X <= x)",
     col=plot.col["CDF"],
     lty=line.types["CDF"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)
par(par.ori)
par(mfrow=c(1,1))
```


## Distribution binomiale en $j$

Quand la probabilité de succès est proche de 1, la distribution et monotone croissante et on la qualifie de ***distribution en forme de $j$***. 

```{r binomial_distrib_j-shaped, fig.width=9, fig.height=4, echo=FALSE, fig.cap="Distribution binomiale en forme de j.  "}

## Drawing parameters
par(mfrow=c(1,2))
par(las=1)

p <- 0.98
binomial.mean <- p*n

xmax <- 5 + 4*binomial.mean

plot(x.shifted,
     dbinom(x,n,p),
     main=paste("Binomial probability (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X=x)",
     col=plot.col["probability"],
     lty=line.types["probability"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)

plot(x.shifted,
     pbinom(x,n,p),
     main=paste("CDF (p=",format(p),", n=", format(n), ")"),
     ylim=c(0,1),
     xlab="X (successes)",
     ylab="P(X <= x)",
     col=plot.col["CDF"],
     lty=line.types["CDF"],
     type="s",
     lwd=2, las=1,
     panel.first=c(abline(v=x.grid,col=plot.col['grid'],lty=line.types["grid"]),
                   abline(h=y.grid,col=plot.col['grid'],lty=line.types["grid"]))
)
par(par.ori)
par(mfrow=c(1,1))
```



## Exemples d'applications de la binomiale

1. **Jeu de dés**: nombre de 6 observés lors d'une série de 10 tirages. 
2. **Alignement de séquences**: nombre d'identités entre deux séquences alignées sans gap.
3. **Analyse de motifs**: nombre d'occurrences d'un motif dans un génome.

**Note: ** le recours à la binomiale présuppose un schema de Bernoulli. Pour les exemples 2 et 3 ceci revient à considérer que les nucléotides se succèdent de façon indépendante, ce qui est assez peu réaliste. 

## Loi de Poisson

La loi de Poisson décrit la probabilité du nombre d'occurrences d'un événement pendant un intervalle de temps fixé, en supposant que le nombre moyen d'événements par unité de temps est constant, et que les événements sont indépendants (les réalisations précédentes n'affectement pas la probabilité des occurrences suivantes).

### Fonction de masse

$$P(X = x) = \frac{\lambda^x}{x!}e^{-\lambda}$$

- $x$ est le nombre d'événements observés;
- $\lambda$ (lettre grècque "lambda") représente l'espérance, autrement dit la moyenne attendue pour le nombre d'événements;
- $e$ est la base de l'exponentielle ($e = `r round(digits=3, exp(1))`$).

## Propriétés de la distribution de Poisson

- **Espérance** (nombre de succès attendus au hasard): $<X> = \lambda$ (par construction)

- **variance**: $\sigma^2 = lambda$ (***variance égale à la moyenne!***)
    
- **Ecart-type**: $\sigma = \sqrt{\lambda}$



## Application : mutagenèse

- On soumet une population à un mutagène (agent chimique, irratiations). Chaque individu subit un certain nombre  de mutations. 

- En tenant compte de la dose de mutagène (temps d'exposition, intensité/concentration), on peut estimer empiriquement le nombre moyen de mutations par individu (*espérance*, $\lambda$).

- La loi de Poisson peut être utilisée pour décrire la probabilité d'observer un nombre donné de mutations ($x=0, 1, 2, ...$).

## Expérience historique de Luria-Delbruck (1943)

En 1943, Salvador Luria et Max Delbruck démontrent que quand on soumet des bactéries en culture à un antibiotique, les mutations qui confèrent la résistance ne sont pas induites par l'antibiotique, mais pré-existent. Leur démonstration repose sur le fait que le nombre de cellules résistantes suit une loi de Poisson ([Luria & Delbruck, 1943, Genetics 28:491–511](https://www.ncbi.nlm.nih.gov/pubmed/17247100)).


## Convergence de la loi binomiale vers la Poisson

<font color="red">A FAIRE</font>


## Exercices

- [html](04_distributions_discretes_exercices.html)
- [pdf](04_distributions_discretes_exercices.pdf)
- [Rmd](04_distributions_discretes_exercices.Rmd)
