Chapter 3:
Probability Concepts & Distributions
Probability and randomness
Probability and randomness are placeholders for incomplete knowledge.
After I shuffled a deck of cards, you might consider the identity of the top card to be “random”.
But is it really?
If you knew the starting positions of the cards and a good HD video of my shuffling, you could surely know the positions of the cards, and which is on top.
Likewise for rolling a die. If we know everything about the starting position, how it was thrown, the texture of the surface, humidity, etc., could we predict what it would roll?
Blood Group Example
Two systems for categorising blood are:
- the Rh system (Rh+ and Rh–)
- the Kell system (K+ and K–)
For any person, their blood type in any one system
is independent of their blood type in any other.
For Europeans in New Zealand,
about 81% are Rh+ and about 8% are K+.
From the table:
If a European New Zealander is chosen at random, what is the probability that they are (Rh+ and K+) or (Rh– and K–)?
- 0.0648 + 0.1748 = 0.2396
Suppose that a murder victim has a bloodstain on him with type (Rh– and K+), presumably from the assailant. What is the probability that a randomly selected person matches this type?
- 0.0152
Discrete probability distributions
Consider the number of eggs \((X)\) in an Adelie penguin’s nest. The values range from \(1\) to \(5\), each with a certain probability (or relative frequency) of occurrence.
- Note the probabilities add to \(1\) because \({1,2,3,4,5}\) is a complete sample space.
The population mean \(\mu_X\) is simply the sum of each outcome multiplied by its probability.
\[\mu_X = E(X)= \sum xP(X=x)=\sum xP(x)\]
In R,
The population variance is given by
\[Var(X)= \sigma_X^2=\sum (x-\mu_X)^2 P(x)\]
The population SD is simply the square-root of the variance.
In R,
Binomial distribution
Consider a variable that has two possible outcomes
(say success and failure, with 50% probabilty each).
This can be described as a “Bernoulli” random variable.
A “Binomial” is just a collection of Bernoulli trials.
Let \(X\) be the number of heads when two coins are tossed.
The count of the number of successes \(X\) out of a fixed total of
\(n\) independent trials follows the binomial distribution.
That is, \(X \sim Bin(n, p)\), where \(p\) the probability of a success.
The binomial probability function \(P(X=x)\) or \(P(x)\)
is given by \[P(x)={n \choose x}p^{x}(1-p)^{n-x}\]
For \(n=10\), \(p=0.3\), the binomial probabilities,
\(P(x)\) for \(x=0,1,2, \dots, 10\), are plotted to the right.
If each of 10 basketball shots succeeded with probability 0.3, this describes the probability of your total score out of 10.
Poisson example
Consider the number of changes that accumulate along a
stretch of a neutrally evolving gene over a given period of time.
This is a Poisson random variable with a
population mean of \(\lambda=kt\), where
\(k\) is the number of mutations per generation, and
\(t\) is the time in generations that has elapsed.
Assume that \(k = 1\times10^{-4}\) and \(t = 500\).
For \(\lambda=kt=0.05\), the Poisson probabilities are shown in the following plot.
What is the probability that at least one mutation has occurred over this period?
\(P(x > 0)=1-P(x = 0)\) is found in R as follows:
Continuous probability distributions
For example, consider a random proportion \((X)\) between \(0\) and \(1\). Here \(X\) follows a (standard) continuous uniform distribution whose (probability) density function \(f(x)\) is defined as follows:
\[f(x)=\begin{cases}{1}~~~\mathrm {for} \ 0\leq x\leq 1,\\[9pt]0~~~\mathrm {for} \ x<0\ \mathrm {or} \ x>1\end{cases}\] This constant density function is the simple one in the graph to the right.
Continuous probability distributions
The density is the relative likelihood of any value of \(x\); that is, the height of the Probability Density Function (PDF). Say, the leaves of a particular tree had mean length 20 cm, SD 2.
The black line is the PDF, or \(f(x)\). The orange area underneath the whole PDF is 1.
Continuous probability distributions
The density is the relative likelihood of any value of \(x\); that is, the height of the Probability Density Function (PDF). Say, the leaves of a particular tree had mean length 20 cm, SD 2.
d <- tibble(x = seq(13, 27, by=0.01),
Density = dnorm(x, 20, 2))
p <- ggplot(d) + aes(x, Density) +
geom_hline(yintercept=0) +
geom_area(colour = 1,
fill = "darkorange",
size = 1.1, alpha = .6)
p +
annotate(geom = "path",
x = c(19.3, 19.3, 13),
y = c(0, rep(dnorm(19.3,20,2),2) ),
arrow = arrow(),
colour = "dodgerblue4", size = 1.1)
The black line is the PDF, or \(f(x)\). The orange area underneath the whole PDF is 1.
The density at 19.3 is \(f(19.3) = 0.1876\)).
Continuous probability distributions
The density is the relative likelihood of any value of \(x\); that is, the height of the Probability Density Function (PDF). Say, the leaves of a particular tree had mean length 20 cm, SD 2.
The black line is the PDF, or \(f(x)\). The orange area underneath the whole PDF is 1.
The area under the curve to the left of the value 19.3 is given by the Cumulative Density Function (CDF), or \(F(x)\). It gives the probability that x < 19.3; \(F(19.3) = 0.3632\).
The Normal (Gaussian) Distribution
The Gaussian or Normal Distribution is parameterised in terms of the mean \(\mu\) and the variance \(\sigma ^{2}\) and its Probability Density Function (PDF) is given by
\[f(x)={\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}\left({\frac {x-\mu }{\sigma }}\right)^{2}}\] A Standard Normal Distribution has mean \(\mu=0\) and standard deviation \(\sigma=1\). It has a simpler PDF:
\[f(z)={\frac {1}{ {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}z^{2}}\] If \(X \sim N(\mu, \sigma)\), you can convert the \(X\) values into \(Z\)-scores by subtracting the mean \(\mu\) and dividing by the standard deviation \(\sigma\).
\[Z={\frac {X-\mu }{\sigma }}\]
We often deal with the standard normal because the symmetric bell shape of the normal distribution remains the same for all \(\mu\) and \(\sigma\).
dfn <- tibble(x=seq(-4,4,length=1000),
`f(x)` = dnorm(x),
`F(x)` = pnorm(x))
p1 <- ggplot(dfn) + aes(x=x,y=`f(x)`) + geom_line() +
geom_vline(xintercept = 0) +
labs(title = "Standard Normal Density",
x = "standard normal deviate, z")
p2 <- ggplot(dfn) + aes(x=x,y=`F(x)`) + geom_line() +
geom_vline(xintercept = 0) +
labs(title = "Cumulative Standard Normal Density",
x = "standard normal deviate, z")
p1/p2 
Example of a normal
The weight of an individual of Amphibola crenata, a marine snail,
is normally distributed with a mean of \(40g\) and variance of \(20g^2\).
What is the probability of getting a snail that weighs between \(35g\) and \(50g\)?
What is the probability of getting a snail that weighs below \(35g\) or over \(50g\)?
Areas (probabilities) under the standard normal
Under standard normal, the areas under the PDF curve are shown below for various situations.
Log-normal distribution
A random variable \(X\) is log-normally distributed, when \(log_e(X)\) follows normal.
Alternatively, if \(X\) follows normal, then \(e^X\) follows log-normal.
R function dlnorm() gives the log-normal density.
Mean & variance:
\[ \begin{align} \mu_X&=e^{\left(\mu +{\frac {\sigma ^{2}}{2}}\right)} \\ \sigma_X^2&=(e^{\sigma ^{2}}-1) e^{(2\mu +\sigma ^{2})} \end{align} \]
dfln <- tibble(x=seq(0, 4, length=1000),
`f(x)` = dlnorm(x),
`F(x)` = plnorm(x))
p1 <- ggplot(dfln) + aes(x=x,y=`f(x)`) + geom_area(alpha=.8) +
labs(title = "Standard log-normal density",
x = "standard log-normal deviate, z")
p2 <- ggplot(dfln) + aes(x=x,y=`F(x)`) + geom_line() +
labs(title = "Cumulative standard log-normal density",
x = "standard log-normal deviate, z")
p1/p2 
Weibull distribution
The PDF of the Weibull distribution is:
\[f(t;\eta,\beta) = \begin{cases} \frac{\beta}{\eta}\left(\frac{x}{\eta}\right)^{\beta-1}e^{-(x/\eta)^{\beta}} ~~x\geq0 ,\\ 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ x<0, \end{cases}\]
\(\beta~(> 0)\) is the called the shape parameter and \(\eta~(> 0)\) is the called scale parameter.
The Weibull distribution becomes the exponential distribution for \(\beta=1\).
The scale parameter \(\eta\) is called the characteristic life because \(\eta\) becomes the quantile with slightly less than two-thirds of the population (63%) below it irrespective of the shape \(\beta\).
Gamma distribution
The probability function of the gamma distribution with shape parameter \(\alpha\) and scale parameter \(\beta\) is given below:
\[\displaystyle {\begin{aligned}f(x)={\frac {\beta ^{\alpha }x^{\alpha -1}e^{-\beta x}}{\Gamma (\alpha )}}\quad {\text{ for }}x>0\quad \alpha ,\beta >0,\\[6pt]\end{aligned}}\] where \(\displaystyle \Gamma (\alpha)=\int _{0}^{\infty }x^{\alpha-1}e^{-x}\,dx.\)
Beta distribution
The beta distribution is bounded on the interval \([0, 1]\) and parameterised by two positive shape parameters, say \(\alpha\) and \(\beta\).
Probability function of the beta distribution:
\[\begin{aligned}f(x;\alpha ,\beta ) ={\frac {x^{\alpha -1}(1-x)^{\beta -1}}{\displaystyle \int _{0}^{1}u^{\alpha -1}(1-u)^{\beta -1}\,du}}={\frac {1}{\mathrm {B} (\alpha ,\beta )}}x^{\alpha -1}(1-x)^{\beta -1}\end{aligned} \] where \(\mathrm {B} (\alpha ,\beta )=\frac {\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}\). Here’s a plot of beta density for various shape parameter combinations.
When \(\alpha=\beta=1\), the beta distribution becomes the continuous uniform distribution.
Small sample effect
For small samples, the shape might be difficult to judge.
set.seed(1234)
dfm <- data.frame(
x=rnorm(50,
mean=80,
sd=12)
)
p1 <- ggplot(dfm) +
geom_histogram(
aes(x=x, y=after_stat(density)),
colour=1
) +
stat_function(
fun = dnorm,
args = list(mean = 80, sd = 12),
geom = "line"
) +
xlim(min(dfm), max(dfm))
p2 <- ggplot(dfm) + aes(x) +
geom_boxplot() +
xlim(min(dfm), max(dfm)) +
theme_void()
library(patchwork)
p1 / p2 + plot_layout(heights = c(5, 1))
Normal quantile plots
In a normal quantile plot, the quantiles of the sample are plotted against the theoretical quantiles of the fitted normal distribution.
The points should roughly lie on a straight line
We can also compare the empirical and theoretical CDFs.
TV viewing time data
Skewed data
The number of people who made use of a recreational facility is in the rangitikei dataset.
The deviation of points from the line indicate that this variable does not conform to a normal distribution.
Fitting lognormal
R package fitdistrplus can be used to estimate the lognormal parameters for people data.
The likelihood function is based on the joint probability of the observed data as a function of the distributional parameters.
The ML method maximises the likelihood function to estimate the distributional (model) parameters. Examine the fit using graphical displays.
Gamma fit
Gamma distribution also fits OK but the outlier remains.
Outliers and subgroups will affect the Maximum Likelihood (ML) employed in the fitdistrplus package.
Fitting of the distribution ' gamma ' by maximum likelihood
Parameters:
estimate Std. Error
shape 1.14299006 0.249453814
rate 0.01593588 0.004314727
Cullen and Frey plot
This is a plot of kurtosis, a measure of excessive peakedness, against the square of a measure of skew.
summary statistics
------
min: 4 max: 470
median: 46
mean: 71.72727
estimated sd: 86.28089
estimated skewness: 3.3264
estimated kurtosis: 17.11618 
Symmetrical distributions such as normal (read corresponding to zero skew) are poor fits. Lognormal also does not fare well.
When the observed data is a mixture from two or more distributions or contain a large number of unusual observations, no single distributional fit will be satisfactory.
TV watching time data
Normal fit is supported.
summary statistics
------
min: 490 max: 2799
median: 1760.5
mean: 1729.283
estimated sd: 567.913
estimated skewness: -0.004088284
estimated kurtosis: 2.32614 
Which distribution is most appropriate?
Remember, theoretical distributions aren’t real—they’re just models—but they can be useful. Keep your purpose in mind.
Choose the simplest distribution that provides an adequate fit.
Data may be best served by a mixture of two or more distributions rather than a single distribution.
