--- title: "Discrete distributions" author: "Jacques van Helden" date: "`r Sys.Date()`" output: beamer_presentation: colortheme: dolphin fig_caption: yes fig_height: 6 fig_width: 7 fonttheme: structurebold highlight: tango incremental: no keep_tex: no slide_level: 2 theme: Montpellier toc: yes revealjs::revealjs_presentation: theme: night transition: none self_contained: true css: ../slides.css html_document: fig_caption: yes fig_height: 6 fig_width: 7 highlight: tango incremental: no keep_md: yes smaller: yes theme: cerulean toc: yes toc_depth: 3 toc_float: yes widescreen: yes ioslides_presentation: css: slides.css fig_caption: yes fig_height: 6 fig_width: 7 highlight: tango smaller: yes toc: yes widescreen: yes pdf_document: fig_caption: yes highlight: zenburn toc: yes toc_depth: 3 slidy_presentation: fig_caption: yes fig_height: 6 fig_width: 7 highlight: tango incremental: no keep_md: yes smaller: yes theme: cerulean toc: yes widescreen: yes font-import: http://fonts.googleapis.com/css?family=Risque subtitle: Probabilities and statistics for biology (CMB STAT1 - STAT2) font-family: Garamond transition: linear --- ```{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) ``` ## Discrete distributions of probabilities The expression ***discrete distribution*** denotes probability distribution of variables that only take discrete values (by opposition to continuous distributions). **Notes:** - In probabilities, the observed variable ($x$) usually represents the number of successes of a series of tests, or the counts of some observation. In such cases, its values are natural numbers ($x \in \mathbb{N}$). - The probability $\operatorname{P}(x)$ takes real values comprised between 0 and 1, but its distribution is said *discrete¨since it is only defined fora set of discrete values of $X$. It is generally represented by a step function. ## Geometric distribution **Application:** waiting time until the first appeearance of an event in a Bernoulli schema. **Examples:** - In a series of dices rollings, count the number rolls ($x$) before the first occurrence of a 6 (this occurrence itself is not taken into account). - Length of a DNA sequence before the first occurrence of a cytosine. ## Mass function of the geometric distribution The ***Probability Mass Function*** (***PMF***) indicates the probability to observe a particular result. For the geometric distribution, it indicates the probability to observe exactly $x$ failures before the first success, in a series of independent trials with a probability of success $p$. $$\operatorname{P}(X = x) = (1-p)^x \cdot p$$ Justification: - The probability of failure for the first trial is $q = 1-p$ (complementary events). - Bernoulli schema $\rightarrow$ the trials are independent $\rightarrow$ the probability of the series is the product of probabilities of its successive outcomes. - One thus computes the product of probabilities of the $x$ initial failures and of the success at the $(x+1)^\text{th}$ trial. *Note: the PMF of discrete distributions relates to the concept of **density** used for continuous distributions.* ## Geometric PMF ```{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)) ``` ## Distribution tails and cumulative distribution function The ***tails*** of a distribution are the areas comprised under the density curve up to a given value (***left tail***) or staring from a given value (***right tail***). - The ***right tail*** indicates the probability to observe a result ($X$) **smaller than or equal to** a given value ($x$): $\operatorname{P}(X \le x$). - **Definition: ** the ***Cumulative Density Function*** (***CDF***) $\operatorname{P}(X \le x)$ indicates the probability for a random variable $X$ to take a value smaller than or equal to a given value ($x$). It corresponds to the left tail of the distribution (including the $x$ value). - The **left tail** of a distribution indicates the probability to observe a result **higher than or equal to** a given value: $\operatorname{P}(X \ge x$). - Note: in the next chapters we will see the use of the right tail of different distributions to measure the ***P value***, in the context of functional enrichment tests, motif over-representations, peak detection, $\ldots$ ## Distribution tails and cumulative distribution function ```{r goemetric_CDF, fig.width=7, fig.height=5.5, out.width="60%", echo=FALSE, fig.cap="**Tails and Cumulative Density Function of the geometric distribution**. "} 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)) ``` ## Binomial distribution The **binomial distribution** indicates the probability to observe a given number of successes ($x$) in a series of $n$ independent trials with constant success probability $p$ (Bernoulli schema). **Binomial PMF** $$\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}$$ **Binomial CDF** $$\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}}$$ **Properties** - **Expectation** (number of successes expected by chance): $ = n \cdot p$ - **Variance**: $\sigma^2 = n \cdot p \cdot (1 - p)$. - Note: the variance of the binomial is inferior to its mean - **Standard deviation**: $\sigma = \sqrt{n \cdot p \cdot (1 - p)}$ ## $i$-shaped binomial distribution The binomial distribution can take various shapes depending on the values of its parameters (success probability $p$, and number of trials $n$). When the expectation ($p \cdot n$) is very small, the binomial distribution is monotonously decreasing and is qualified of ***$i$-shaped***. ```{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)) ``` ## Asymmetric bell-shaped binomial distribution When the probability is relatively high but still lower than $0.5$, the distribution takes the shape of an asymmetric bell. ```{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)) ``` ## Symmetric bell-shaped binomial When the success probability $p$ is exactly $0.5$, the binomial distribution takes the shape of a symmetrical bell. ```{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)) ``` ## $j$-shaped binomial distribution Then the success probability is close to 1, the distirbution is monotonously increasing and is qualified of ***$j$-shaped distribution. ```{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)) ``` ## Examples of applications of the binomial 1. **Dices**: number of 6 observed during a series of 10 dice rolls 2. **Sequence alignment**: number of identities between two sequences alignmed without gap and with an arhbitrary offset. 3. **Motif analysis**: number of occurrences of a given motif in a genome. **Note: ** the binomial assumes a Bernoulli schema. Forexamples 2 and 3 this amounts to consider that nucleotides are concatenated in an independent way, which is quite unrealistic. ## Poisson law The Poisson law describes the probability of the number of realisations of an event during a fixed time interval, assuming that the average number of events is constant, and that the events are independent (previous realisations do not affect the probabilities of future realisations). ### Poisson Probability Mass Function $$P(X = x) = \frac{\lambda^x}{x!}e^{-\lambda}$$ - $x$ is the number of event realisations - $\lambda$ (Greek letter "lambda") represents the expectation, i.e. the average number of occurrences that would be obtained by running the same test an infinite number of times; - $e$ is the exponential base ($e = `r round(digits=3, exp(1))`$). ## Properties of the Poisson distribution - **Expectation** (number of realisations expected by chnace): $ = \lambda$ (by construction) - **variance**: $\sigma^2 = lambda$ (***the variance equals the mean!***) - **Standard deviation**: $\sigma = \sqrt{\lambda}$ ## Application: mutagenesis - A bacterial population is submitted to a mutagen (chemical agent, irradiations). Each cell is affected by a particular number of mutations. - Taking into account the dosis of the mutagen (exposure time, intensity, concentration) one could take an empirical measure of the mean number of mutations by individual (expectation, $\lambda$). - The Poisson law can be used to describe the probability for a given cell to have a given number of mutations ($x=0, 1, 2, ...$). ### Historical experiment by Luria-Delbruck (1943) In 1943, Salvador Luria and Max Delbruck demonstrated that when cultured bacteria are treated by an antibiotic, the mutations that confer resistance are not induced by the antibiotic itself, but preexist. Their demonstration relies on the fact that the number of antibiotic-resistant cells follows a Poisson law ([Luria & Delbruck, 1943, Genetics 28:491–511](https://www.ncbi.nlm.nih.gov/pubmed/17247100)). ## Convergence of the binomial towards the Poisson Under some circumstances, the binmial law converges towards a Poisson. - very small probability of success ($p \ll 1$) - large number of trials ($n$) TO DO ## Netative binomial: number of successes before the $r^{th}$ failure The **negative binomial** distribution (also called **Pascal distribution**) indicates the probability of the number of successes ($k$) before the $r^{th}$ failure, in a Bernoulli schema with success probability $p$. $$\mathcal{NB}(k|r, p) = \binom{k+r-1}{k}p^k(1-p)^r$$ This formula is a simple adaptation of the binomial, with the difference that we know that the last trial must be a failure. The binomial coefficient is thus reduced to choose the $k$ successes among the $n-1 = k+r-1$ trials preceding the $r^{th}$ failure. ## Negative binomial: alternative formulations It can also be adapted to indicate related probabilities. - Number of **failures** ($r$) before the $k^{th}$ **success**. $$\mathcal{NB}(r|k, p) = \binom{k+r-1}{r}p^k(1-p)^r$$ - Number of **trials** ($n=k+r-1$) before the $r^{th}$ **failure**. $$\mathcal{NB}(n|r, p) = \binom{n-1}{r-1}p^{n-r}(1-p)^r$$ ## Negative binomial density ```{r negbin, echo=FALSE, fig.width=8, fig.height=4, fig.cap="Negative binomial. "} k <- 0:25 # Number of successes par(mfrow=c(2,2)) par(mar=c(4,4,1,1)) for (r in c(1, 2, 5, 10)) { P.k <- dnbinom(k, size = r, prob = 0.75) plot(k, P.k, type="h", panel.first = grid(), lwd=3, ylim=c(0,1), las=1, xlab="Number of successes (k)", ylab="P(K=k)", col="darkblue") legend("topright", legend=paste("Last failure: r = ", r)) } par(mfrow=c(1,1)) ``` ## Properties of the negative binomial The variance of the negative binomial is higher than its mean. It is therefore sometimes used to model distributions that are over-dispersed by comparisong with a Poisson. $$\mathcal{NB}(r|k, p) = \binom{k+r-1}{r}p^k(1-p)^r$$ - Parameters: - $p$: probability of success at each trial - $r$: number of failures - $k$: number of successes before the $r^{th}$ failure - Mean: $\frac{pr}{1-p}$ - Variance: $\frac{p(1-p)}{p^2}$ ## Exercise -- Negative binomial Each student chooses a value for the maximal number of failures ($r$). 1. Read carefully the help of the negative binomial functions: `help(NegBinomial)` 2. **Random sampling**: draw of $rep=100000$ random numbers from a negative binomial distribution (`rndbinom()`) to compute the distribution of the number of successes ($k$) before the $r^{th}$ failure. 3. Compute the expected mean and variance of the negative binomial. 4. Compute the mean and variance from your sampling distribution. 5. Draw an histogram with the number of successes before the $r^{th}$ failure. 6. Fill up the form on the [collective result table](https://docs.google.com/spreadsheets/d/1Kl_0ln0_dZycK17Nqyu44kw9R0dtVp5lflXRtN7pAhA/edit#gid=0) ## Solution to the exercise -- negative binomial ```{r out.width="50%"} r <- 6 # Number of failures p <- 0.75 # Failure probability rep <- 100000 k <- rnbinom(n = rep, size = r, prob = p) max.k <- max(k) exp.mean <- r*(1 - p)/p rand.mean <- mean(k) exp.var <- r*(1 - p)/p^2 rand.var <- var(k) hist(k, breaks = -0.5:(max.k + 0.5), col = "grey", xlab = "Number of successes (k)", las = 1, ylab = "", main = "Random sampling from negative binomial") abline(v = rand.mean, col = "darkgreen", lwd = 2) abline(v = exp.mean, col = "green", lty = "dashed") arrows(rand.mean, rep/20, rand.mean + sqrt(rand.var), rep/20, angle = 20, length = 0.1, col = "purple", lwd = 2) text(x = rand.mean, y = rep/15, col = "purple", labels = paste("sd =", signif(digits = 2, sqrt(rand.var))), pos = 4) legend("topright", legend = c( paste("r =", r), paste("exp.mean =", signif(digits = 4, exp.mean)), paste("mean =", signif(digits = 4, rand.mean)), paste("exp.var =", signif(digits = 4, exp.var)), paste("var =", signif(digits = 4, rand.var)) )) kable(data.frame(r = r, exp.mean = exp.mean, mean = rand.mean, exp.var = exp.var, var = rand.var), digits = 4) ``` **************************************************************** # Negative binomial for over-dispersed counts ## Exercises - [html](04_distributions_discretes_exercices.html) - [pdf](04_distributions_discretes_exercices.pdf) - [Rmd](04_distributions_discretes_exercices.Rmd)