Week 5 - Univariate Probability Distributions

library(ggplot2)
library(MASS)
library(PASWR2)

Objectives

• Compute probabilities related to discrete and continuous random variables
• Use experimental and theoretical probabilities to formulate and solve applied problems

Motivation:

Discrete distributions can be used to model the number of failures until a successful rocket launch, the number of passing students in a class, or the number of taxis that pass a street corner, as well as many other phenomena with countable outcomes. Continuous distributions are used to model measurement variables such as weight, height, and time.

Motivational Problem (Pr. 20 / page 309 in text)

The weekly production of a banana plantation can be modeled with a normal random variable that has a mean of 5 tons and a standard deviation of 2 tons.

1. Find the probability that, in at most 1 out of the 8 randomly chosen weeks, the production has been less than 3 tons.
2. Find the probability that at least 3 weeks are needed to obtain a production greater than 10 tons.

• Solution towards the end of the notes!

Discrete Univariate Distributions

4.2.1 Discrete Uniform Distribution

The random variable X is said to follow a discrete uniform distribution with parameter n (where \(n ∈ ℕ)\) if the probability X takes on the value x is the same for all x, where \(x = x_1, x_2, …, x_n\).

For Discrete Uniform Distribution one has:

1. \[P(X=x_i|n)= \frac{1}{n}\] for \(i=1,2,...,n\)

2. \[E[X]=\frac{1}{n} \sum_{i=1}^n x_i\]

3. \[Var[X]=\frac{1}{n} \sum_{i=1}^n (x_i-E[X])^2\]

Exercise: When \(x_i = i\) for \(i = 1, …, n\), it can be shown that \(E[X] = \frac{n+1}{2}\) and that \(Var[X] = \frac{n^2−1}{12}\), respectively.

Example 4.1

One light bulb is randomly selected from a box that contains a 40-watt light bulb, a 60-watt light bulb, a 75-watt light bulb, a 100-watt light bulb, and a 120-watt light bulb. Write the probability function for the random variable that represents the wattage of the randomly selected light bulb, and determine the mean and variance of that random variable.

Solution: The random variable X can assume the set of values \(Ω = {40, 60, 75, 100, 120}\). The probability density function for the random variable X is

\(P(X=x|n=5)=\frac{1}{5}\) for \(x=40, 60, 75, 100, 120\)

Try the code

Watts <- c(40, 60, 75, 100, 120)
n <- length(Watts)
meanWatts <- (1/n)*sum(Watts)
varWatts<- (1/n)*sum((Watts - meanWatts)^2)
ans <- c(meanWatts, varWatts)
ans

4.2.2 Bernoulli and Binomial Distributions

A Bernoulli trial is a random experiment with only two possible outcomes. The outcomes are mutually exclusive and exhaustive;e.g., success or failure, true or false, alive or dead, male or female, etc. A Bernoulli random variable, X, can take on two values, where X(success) = 1 and X(failure) = 0. The probability that X is a success is \(\pi\), and the probability that X is a failure is \(1-\pi\).

For Bernoulli distribution pdf, mean, and variance of a Bernoulli random variable are:

\[X \sim Bernoulli(\pi)\]

1. \[P(X=x|\pi)= \pi^x(1-\pi)^{1-x}\] for \(x=0,1\)

2. \[E[X]=\pi\]

3. \[Var[X]=\pi(1-\pi)\]

When a sequence of Bernoulli trials conforms to the following list of requirements, it is called a binomial experiment:

1. The experiment consists of a fixed number (n) of Bernoulli trials.
2. The probability of success for each trial, denoted by \(\pi\), is constant from trial to trial. The probability of failure is \(\rho = 1-\pi\)
3. The trials are independent.
4. The random variable of interest, X, is the number of observed successes during the n trials.

For Binomial distribution pdf, mean, and variance of a binomial random variable are:

\[X \sim Bin(n, \pi) \]

1. \[P(X=x|n, \pi)= {n \choose x} \pi^x(1-\pi)^{n-x}\] for \(x=0,1, ..., n\)

2. \[E[X]=n \pi\]

3. \[Var[X]=n \pi(1-\pi)\]

Exercise (optional): Prove 1,2 and 3 above.

Code below is used to create graphs that represent the probability density function and the cumulative distribution function for a Bin(8, 0.3) random variable:

Try it!

opar <- par(no.readonly = TRUE)
par(mfrow=c(1, 2), pty = "s")
x <- 0:8
px <- dbinom(x = x, size = 8,prob =  0.3)

# plot probability vs corresponding x-values
plot(x, px, type = "h", xlab = "x", ylab="P(X = x)",
    main = "PDF of X~Bin(8, 0.3)")
xs <- rep(0:8, round(dbinom(0:8, 8, .3)*100000, 0)) # creates vector of 0s, 1s, ...8s, according to the expected number of each one out of 10,000 trials

# plot the empirical cdf
plot(ecdf(xs), main = "CDF of X~Bin(8, 0.3)",
    ylab = expression(P(X<=x)), xlab = "x")
par(opar)

Q: Change the probability, \(\pi\) parameter to value close to zero. What do you observe? How does it skew the distribution. Now, increase `n’ to bigger value, even bigger. Do you notice correction of the skewness?

Exercise: Graph the pdf and cdf for \(X \sim Bin(10,\frac{1}{2})\) (# of heads of tossing a fair coin 10 times experiment) it is part of the weekly programming exercises Rmd file.

R functions for Binomial Distribution Calculations

Description

Density, distribution function, quantile function and random generation for the binomial distribution with parameters size and prob.

This is conventionally interpreted as the number of ‘successes’ in size trials.

Usage

dbinom(x, size, prob, log = FALSE) this computes the pdf values

pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE) this gives \(P(X \le x)\) the pdf values

qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE) this gives the quantiles \(P(X \le k) \ge p\)

rbinom(n, size, prob)

Arguments

x, q - vector of quantiles. p - vector of probabilities. n - number of observations. size - number of trials (zero or more). prob - probability of success on each trial.

Example: Using the function rbinom(), one can generate 1000 samples of a \(Bin(n = 5, \pi = 0.5)\) as shown below:

set.seed(123)
x <- rbinom(1000, 5, .5)
table(x)/1000 # The empirical distribution

Example 4.3

Binomial Calculation Consider the problem of calculating the probability of obtaining 6 or more heads in 10 tosses of a weighted coin (not fair), where the probability of obtaining a head in any given trial is 0.33.

# adding probabilities 
sum(dbinom(x = 6:10, size = 10, prob = 0.33))          # P(X >= 6)=P(6)+P(7)+P(8)+P(9)+P(10)

# or
1 - pbinom(5, 10, 0.33)              # 1 - P(X <= 5)
#or 
pbinom(5, 10, 0.33, lower = FALSE)   # P(X >= 6)
# or 
1 - sum(dbinom(5:0, 10, 0.33)) # 1 - P(X <= 5)

4.2.3 Poisson Distribution

The Poisson distribution is very popular for modeling the number of times particular events occur in given times or on defined spaces, or in general unit interval of time (units of time may vary).

E.g.

• number of phone calls to 911 between 1 a.m. and 2 a.m.

• number of accidents at a busy street corner during a 24-hour period

• number of typographical errors on a single page of this book

See definition of Poisson Process in the text, page 256.

Suppose there is an experiment that satisfies the three criteria for an approximate Poisson process with parameter \(\lambda\), \(\lambda\) mean number of occurrences on a unit interval of time.

Let X represent the number of outcomes in an interval of length 1 (t = 1).

Then one has:

Poisson Distribution mean, and variance of a binomial random variable are:

\[X \sim Pois(\lambda)\]

1. \[P(X=x|\lambda)=\frac{\lambda^xe^{-\lambda}}{x!}\] for \(x=0,1,2,...\)

2. \[E[X]=\lambda\]

3. \[Var[X]=\lambda\]

Review the code to see how to represent a pdf and cdf for a Pois(λ = 1) random variable. Although the values x can take on with the Poisson distribution are $0, 1, 2,…, $ the probability a Poisson random variable has a value of eight or greater when λ = 1 is extremely small (< 0). Consequently, the pdf, and cdf do not extend beyond eight.

opar <- par(no.readonly = TRUE)
par(mfrow=c(1, 2), pty = "s")

x <- 0:8
px <- dpois(x,lambda = 1) # computes the pdf at the points x

# plot probability vs x-values
plot(x, px, type = "h", xlab = "x", ylab="P(X = x)",
    main = "PDF of X ~ Pois(1)")

xs <- rep(0:8, round(dpois(0:8, 1)*100000, 0)) # create vector of 100,000 values according to the expected number of occurrences of 0s,1s, 2s, ...and we cut at 8s other values can be ignored
# plot the empirical cdf
plot(ecdf(xs), main = "CDF of X ~ Pois(1)",
    ylab = expression(P(X <=x)), xlab = "x")

par(opar)

Note: The Poisson distribution is particularly appropriate for modeling “rare” phenomena or outcomes where the probability of success is small.

Given n independent Poisson random variables \(X_1, X_2,…, X_n\) with parameters \(λ_1, λ_2,…, λ_n\) respectively, \(Y= \sum_{i=1}^n X_i \sim Pois(\sum_{i=1}^n \lambda_i=\lambda)\)

Example 4.6 More accidents are registered in auto body repair shops during the months of May and June than in the rest of the year. Suppose a particular auto body repair shop has an average of four accidents per month. What is the probability there will be more than seven accidents in this auto body shop during the month of May? What is the probability no more than three accidents will occur during the months of May and June?

Solution: Assuming accidents in the auto body shop follow an approximate Poisson process, the probability of x accidents in one month is

\[P(X=x|4)=\frac{4^xe^{-4}}{x!}\] for \(x=0,1,2,...\)

Using R:

# lower parameter TRUE - computes probabilities up to the point, FALSE to the right of the point
1 - ppois(q = 7, lambda = 4, lower = TRUE)               # P(X > 7|4) = 1 - P(X<=7|4)

ppois(q = 7, lambda = 4, lower = FALSE)    # P(X > 7|4) notice lower = FALSE

ppois(q = 3, lambda = 8) # P(X <= 3|8)

R functions for the Poisson Distribution Calculations

Description

Density, distribution function, quantile function and random generation for the Poisson distribution with parameter lambda.

Usage

dpois(x, lambda, log = FALSE)

ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)

qpois(p, lambda, lower.tail = TRUE, log.p = FALSE)

rpois(n, lambda)

Arguments

x - vector of (non-negative integer) quantiles. q - vector of quantiles. p - vector of probabilities. n - number of random values to return. lambda - vector of (non-negative) means.

4.2.4 Geometric Distribution

A random variable X that counts the number of Bernoulli trials that result in failure before the first success is called a geometric random variable. Clearly, the probability of a success after r failures is \(π × (1 − π)^r\), which leads to the geometric probability distribution function where \(\rho = 1 − π\) is the probability of failure.

The pdf, mean, and variance for a geometric random variable are shown below:

\[X \sim Geo(\pi)\] \[P(X=x|\pi)=\pi\rho^x\] for \(x=0,1,2 ...\)

\[E[X]=\frac{\rho}{\pi}\] \[Var[X]=\frac{\rho}{\pi^2}\]

The code below can be used to create graphs that represent the probability density function and the cumulative distribution function for a \(Geo(π = 0.3)\) random variable.

opar <- par(no.readonly = TRUE)
par(mfrow=c(1, 2), pty = "s")
x <- 0:15
px <- dgeom(x = x, prob = .3)

plot(x, px, type = "h", xlab = "x", ylab="P(X = x)",
     main = "PDF of X ~ Geom(0.3)")
xs <- rep(0:15, round(dgeom(0:15, 0.3)*100000, 0))
plot(ecdf(xs), main = "CDF of X ~ Geom(0.3)",
     ylab = expression(P(X <=x)), xlab = "x")
par(opar)

Pr. 18 / page 308 in text Suppose the percentage of drinks sold from a vending machine are 80% and 20% for soft drinks and bottled water, respectively.

1. What is the probability that on a randomly selected day, the first soft drink is the fourth drink sold?
2. Find the probability that exactly 1 out of 10 drinks sold is a soft drink.

Solution: Let X = number of waters (failures) purchased before the first soft drink is purchased. Then, \(X \sim Geo(0.80)\).

dgeom(x = 3, prob = 0.80) # P(X = 3)

2. Let X = number of soft drinks sold. Then, \(X \sim Bin(10, 0.80)\) and \(P(X = 1) = 0\) since one has:

pbinom(q = 1, size = 10, prob = 0.8) # P(X = 1)

Optional to review: 4.2.5 Negative Binomial Distribution & 4.2.6 Hypergeometric Distribution

4.3 Continuous Univariate Distributions

4.3.1 Uniform Distribution (Continuous)

The continuous uniform distribution has a pdf that is constant over a closed interval.

An important application of the uniform distribution includes random number generation.

Uniform Distribution

\[X \sim Unif(a,b)\]

1. \[f(x|a,b)=\frac{1}{b-a}\] for \(a \le x \le b\)

2. \[E[X]=\frac{b+a}{2}\]

3. \[Var[X]=\frac{(b-a)^2}{12}\]

Generating Pseudo Random Numbers using R command rdist where dist is distribution name. Use runif for the next example.

Example:

For \(X \sim Unif(2,6)\) one has \(E[X]=4\) and \(Var[X]=\frac{16}{12}=\frac{4}{3}\)

Estimate \(E[X]\) and \(Var[X]\) by using a large number of values drawn at random from a Unif(2, 6).

set.seed(13)
TM <- 4 # true mean
TV <- 4/3 # true variance

x <- runif(n = 10000,min = 2,max = 6)
Ex <- mean(x)
Ex

Vx <- mean((x-Ex)^2)
Vx

PE <- c(PercentErrorMean  = abs(Ex - TM)/TM*100, 
        PercentErrorVariance = abs(Vx - TV)/TV*100) # percent error of the estimate from the true mean and var
ANS <- c(SimMean = Ex, TheMean = TM, SimVar = Vx, TheVar = TV)
ANS

4.3.2 Exponential Distribution

If W is the waiting time until the first outcome of a Poisson process with mean \(\lambda > 0\) (# occurrences over unit interval), then the pdf, mean, and variance for W are shown below:

1. \[f(x| \lambda)=\lambda e^{-\lambda x}\] for \(x \ge 0\), and 0 otherwise

2. \[E[X]=\frac{1}{\lambda}\]

3. \[Var[X]=\frac{1}{\lambda^2}\]

Consequently, when \(w > 0\), the pdf of W is \(F'(w) = f(w) = λe^{−λw}\).

The exponential distribution is characterized by a lack of memory property and is often used to model lifetimes of electronic components as well as waiting times for Poisson processes.

A random variable is said to be memoryless if

\[P(X > t_2 + t_1 |X > t_1) = P(X > t_2)\] for all \(t_1, t_2\)

Exercise: Try to show that above is equivalent to

\[P(X > t_2 + t_1)=P(X > t_1) P(X > t_2)\] for all \(t_1, t_2\)

R functions for the Exponential Distribution

Description

Density, distribution function, quantile function and random generation for the exponential distribution with rate rate (i.e., mean 1/rate).

Usage

dexp(x, rate = 1, log = FALSE)

pexp(q, rate = 1, lower.tail = TRUE, log.p = FALSE)

qexp(p, rate = 1, lower.tail = TRUE, log.p = FALSE)

rexp(n, rate = 1)

Arguments

x, q - vector of quantiles. p - vector of probabilities. n- number of observations. If length(n) > 1, the length is taken to be the number required. rate - vector of rates. log, log.p - logical; if TRUE, probabilities p are given as log(p). lower.tail - logical; if TRUE (default), probabilities are P[X ≤ x], otherwise, P[X > x].

Example 4.17, page 279 in text

Exponential Distribution: Light Bulbs If the life of a certain type of light bulb has an exponential distribution with a mean of 8 months, find

1. The probability that a randomly selected light bulb lasts between 3 and 12 months.
2. The 95th percentile of the distribution.
3. The probability that a light bulb that has lasted for 10 months will last more than 25 months.

Solution:

YOUR CODE HERE:

# uncomment the code below to see solution - explain why
pexp(q = 12, rate = 1/8) - pexp(q = 3, rate = 1/8)

f1 <- function(x){(1/8) * exp(-x/8)}        # define f1
(integrate(f1, lower = 3, upper = 12)$value)   # integrate f1

# the 95 quantile/percentile
qexp(0.95, 1/8) 

pexp(25, 1/8, lower = FALSE)/pexp(10, 1/8, lower = FALSE)
# or
1 - pexp(15, 1/8)
# or
pexp(15, 1/8, lower = FALSE)

4.3.7 Normal (Gaussian) Distribution

most important distribution in statistical applications

many numerical populations have distributions that can be approximated with the normal distribution.

The pdf, mean, and variance for a normal random variable X with mean μ and variance \(\sigma^2\) are provided: \[X \sim N(\mu, \sigma^2)\] \[ f(x|\mu,\sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}\] for \(-\infty < x <\infty\)

\[E[X]=\mu\]

\[Var[X]=\sigma^2\]

R functions for the Normal Distribution

Description

Density, distribution function, quantile function and random generation for the normal distribution with mean equal to mean and standard deviation equal to sd.

Usage

dnorm(x, mean = 0, sd = 1, log = FALSE)

pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)

qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)

rnorm(n, mean = 0, sd = 1)

Arguments

x, q - vector of quantiles. p - vector of probabilities. n - number of observations. If length(n) > 1, the length is taken to be the number required. mean - vector of means. sd - vector of standard deviations.

Example 4.25

Scores on a particular standardized test follow a normal distribution with a mean of 100 and standard deviation of 10.

1. What is the probability that a randomly selected individual will score between 90 and 115?
2. What score does one need to be in the top 10%?
3. Find the constant c such that \(P(105 \le X \le c) = 0.10\).

# (a)
pnorm(q = 115,mean =  100, sd = 10) - pnorm(q = 90, mean = 100, sd = 10) # see why
# (b) 
qnorm(p = .90, mean = 100, sd = 10) # check why
# (c) 
qnorm(0.10 + pnorm(105, 100, 10), 100, 10) # see the text page 299.

Motivational Problem (Pr. 20 / page 309 in text)

The weekly production of a banana plantation can be modeled with a normal random variable that has a mean of 5 tons and a standard deviation of 2 tons.

1. Find the probability that, in at most 1 out of the 8 randomly chosen weeks, the production has been less than 3 tons.
2. Find the probability that at least 3 weeks are needed to obtain a production greater than 10 tons.

Solution:

1. Let X = number of weeks where production is less than 3 tons. Then, \(X \sim Bin(8, 0.1587)\). Let W = production, \(W \sim N(5, 2)\) and \(P(W \le 3) = 0.1587\). So, P(X ) = 0.6298. Check the code below!

p <- pnorm(3, 5, 2) # P(W < 3)
pbinom(1, 8, p)

2. Let \(X_n\) = number of weeks where production is less than 10 tons before the first week with production over 10 tons. \(X_n \sim Geo(\pi = P(W \ge 10) = 0.0062)\), and \(P(X_n \ge 2) = 1 − P(Xn \le 1) = 0.9876\).

PI <- 1 - pnorm(q = 10,mean =  5, sd = 2) # P(W >= 10)
PI

1 - pgeom(1, PI)
0.9876192

Quantile-Quantile Plots for Normal Distributions (Q-Q Plots)

Many of the techniques and testing assume the underlying distribution is normal.

One of the more useful graphical procedures for assessing distributions is the quantile-quantile plot. (Recall from Section 2.7.3 that this graph is also called a Q-Q plot.)

R function qqnorm() helps determine whether the underlying distribution is normal.

** Example** Consider the values stored in the variable scores of the data frame SCORE (in the PASWR2 package) and shown below which are the scores a random sample of 20 college freshmen received on a standardized test.

119 107 96 107 97 103 94 106 87 112 99 99 90 106 110 99 105 100 100 94

To compute the pairs of values plotted in an R quantile-quantile plot for the variable scores of the data frame SCORE inspect the code below:

n <- length(SCORE$scores) # number of observations

X <- (1:n - 1/2)/n # 100% split into n=20 equal intervals
Xs <- qnorm(X) # quantiles of X
Ys <- sort(SCORE$scores) # sorted scores
plot(Xs, Ys) # plot Ys vs Xs
quantile(Xs, c(0.25, 0.75)) 
quantile(Ys, c(0.25, 0.75))

qqnorm(SCORE$scores) #using qqnorm
qqline(SCORE$scores) #adding line through 1st and 3rd quantile

You can do this using ggplot2 much nicer:

ggplot(data = SCORE, aes(sample = scores)) + stat_qq()

REMINDER: DO NOT MISS TO REVIEW THE PROGRAMMING EXERCISES before you attempt the quiz for the week!