knitr::opts_chunk$set(dev = "ragg_png")
options(digits=5) # set number of digits equal to 5

Packages and functions

# tidyverse
library(tidyverse)

# others
library(gridExtra)
library(kableExtra)
library(glue)
library(highcharter)
# library(scales)

Define a function for printing data frames with kable:

mykable <- function(data) {
  knitr::kable(data) %>% 
    kable_styling()
}

# Colors
colors <- c('darkred', 'deepskyblue', 'darkgreen', 'darkgoldenrod', 'darkblue')

Exercise 1

x <- c(15.58, 15.9, 16, 16.1, 16.2)
p1 <- c(0.15, 0.21, 0.35, 0.15, 0.14)
p2 <- c(0.14, 0.05, 0.64, 0.08, 0.09)

exp_value_1 <- sum(x*p1)
exp_value_2 <- sum(x*p2)

var_1 <- sum(p1*(x-exp_value_1)**2)
var_2 <- sum(p2*(x-exp_value_2)**2)

cat(
  'E[X] method 1:', exp_value_1, 
  '\nE[X] method 2:', exp_value_2, 
  '\nVar(X) method 1:', var_1,
  '\nVar(X) method 2:', var_2
  )

## E[X] method 1: 15.959 
## E[X] method 2: 15.962 
## Var(X) method 1: 0.033979 
## Var(X) method 2: 0.028167

# plot the distribution

pdf_1 <- ggplot() + 
  geom_col(aes(x, p1), fill='gold', colour="black") +
  labs(title="Probability distribution 1", 
       x='x',
       y='p(x)') +
  ylim(0,1)+
  scale_x_continuous(breaks=x, labels=x)+
  theme_bw() 

pdf_2 <- ggplot() + 
  geom_col(aes(x, p2), fill='coral', colour="black") +
  labs(title="Probability distribution 2", 
       x='x',
       y='p(x)') +
  ylim(0,1)+
  scale_x_continuous(breaks=x, labels=x)+
  theme_bw() 


grid.arrange(pdf_1, pdf_2, nrow=1)

Exercise 2

The waiting time, in minutes, at the doctor’s is about 30 minutes, and the distribution follows an exponential pdf with rate 1/30.

A) Simulate the waiting time for 50 people at the doctor’s office and plot the relative histogram

set.seed(1) 
lambda = 1/30
people_wt <- rexp(50, lambda)

ggplot() + 
  geom_histogram(aes(people_wt),  bins = 10, fill='gold', colour="black") +
  labs(title="Simulation of the waiting time for 50 people", 
       x='People Waiting Time (min)',
       y='Count')+
  scale_x_continuous(
    breaks = seq(0, 160, 20),
  ) + 
  theme_bw()

B) What is the probability that a person will wait for less than 10 minutes?

This is given by the cumulative distribution, evaluated at \(t=10\).

cat('Probability that a person will wait for less than 10 minutes:', pexp(10, lambda)*100, '%')

## Probability that a person will wait for less than 10 minutes: 28.347 %

C) Evaluate the average waiting time from the simulated data and compare it with the expected value

From the theory, the expected value of the exponential distribution is \[ E[x] = \frac{1}{\lambda}=30 \]

It is possible to compute it using the definition.

ex <- integrate(function(x){x*dexp(x, lambda)}, lower = -Inf, upper = Inf)$value

cat('Average waiting time by integration:', ex, 'minutes\n')

## Average waiting time by integration: 30 minutes

cat('Error:', abs(ex-30)/30)

## Error: 2.059e-12

Let’s evaluate it from the simulated data.

ex <- mean(people_wt)
cat('Simulated average waiting time:', ex, 'minutes\n')

## Simulated average waiting time: 29.521 minutes

cat('Error:', formatC(abs(ex-30)/30, format = "e"))

## Error: 1.5958e-02

D) What is the probability for waiting more than one hour before being received?

It is \(1-P[T\leq1h]=1-P[T\leq60min]\).

cat('Probability for waiting more than one hour:', (1-pexp(60, lambda))*100, '%')

## Probability for waiting more than one hour: 13.534 %

Exercise 3

Task

Let’s suppose that on a book, on average, there is one typo error every three pages. If the number of errors follows a Poisson distribution, plot the pdf and cdf, and calculate the probability that there is at least one error on a specific page of the book.

Solution

I assume \(\lambda=1/3\), as the expected value of the Poisson distribution is equal to \(\lambda\). The time unit is 1 page.

set.seed(1)

x=0:10
nerrors <- data.frame(X=factor(x, levels = x), pdf=dpois(x, 1/3), cdf=ppois(x, 1/3))

# Plot
pdf <- ggplot() + 
  geom_col(data=nerrors, aes(X, pdf), fill='gold', colour="black") +
  labs(title="Pdf of the number of errors", 
       x='Number of errors',
       y='f(x)') +
  ylim(0,1)+
  theme_bw() 

cdf <- ggplot(nerrors) + 
  geom_point(aes(X, cdf), color='black', size=3) +
  geom_line(aes(X, cdf, group=1), color='black', linetype="dashed") +
  labs(title="Cdf of the number of errors", 
       x='Number of errors',
       y='f(x)') +
  ylim(0,1)+
  theme_bw() 

plot(pdf)

plot(cdf)

The probability that there is at least one error on a specific page of the book is \(1-P[n=0]\), where \(n\) is the number of errors.

cat('Probability that there is at least one error on a specific page of the book:', (1-dpois(0, 1/3))*100, '%')

## Probability that there is at least one error on a specific page of the book: 28.347 %

Exercise 4

Task

We randomly draw cards from a deck of 52 cards, with replacement, until one ace is drawn. Calculate the probability that at least 10 draws are needed.

Solution

The probability of drawing 1 ace from a deck of 52 card is \[ p=\frac{4}{52}=\frac{1}{13} \]

The probability that at least 10 trials are needed until one ace is drawn can be easily achieved using the Binomial distribution. Indeed, it is equal to the probability of having 0 successes in the first 9 trials.

p=1/13

prob_10_a <- dbinom(0, 9, p) 
cat('Probability that at least 10 draws are needed until one ace is drawn:', prob_10_a*100, '%')

## Probability that at least 10 draws are needed until one ace is drawn: 48.657 %

The same result can be achieved in a more complex way using the Negative Binomial distribution.

The probability that exactly \(d\) trials are needed before the first success is given by the Pascal (or Negative Binomial) distribution: \[ P[D=d] = Bneg(r=1|d, p) = p(1-p)^d \] So the probability that at least 10 draws are needed is \[ P[D\geq10] = 1-P[D<10] = 1-\sum_{i=1}^{9}P[D=i] = 1-\sum_{i=1}^{9} Bn(r=1|i, p) \] In R, we can use the function dnbinom(x, size, prob), that represents the number of failures which occur in a sequence of Bernoulli trials before a target number of successes is reached.

’x’ = number of failures before the target number of successes is reached
’size’ = number of successes
’prob’ = p

Using this function, the probability that at least 10 draws are needed can be computed as \[ P[D\geq10] = 1-P[D<10] = 1-\sum_{i=0}^{8} \text{dnbinom}(i, 1, p) \]

prob_10_b <- 1-pnbinom(8, 1, p)

cat('Probability that at least 10 draws are needed until one ace is drawn:', prob_10_b*100, '%')

## Probability that at least 10 draws are needed until one ace is drawn: 48.657 %

Another possible way to solve this exercise is to use the geometric distribution, as the probability we are looking for is the probability of having the first success after 10 or more trials.

prob_10_c <- 1-pgeom(8, p)

cat('Probability that at least 10 draws are needed until one ace is drawn:', prob_10_c*100, '%')

## Probability that at least 10 draws are needed until one ace is drawn: 48.657 %

Exercise 5

The time it takes a student to complete a TOLC-I University orientation and evaluation test follows a density function of the form \[ \begin{equation} f(t) = \begin{cases} c(t-1)(2-t) & 1<t<2\\ 0 & \text{otherwise} \end{cases} \end{equation} \] where \(t\) is the time in hours.

a) Using the integrate() R function, determine the constant \(c\) (and verify it analytically)

To find \(c\), one must impose that the distribution \(f(t)\) is normalized. Analytically: \[ 1 = \int_{-\infty}^{\infty}f(t)dt = \int_{1}^{2}c(t-1)(2-t)dt = c\left[-\frac{t^3}{3}+3\frac{t^2}{2}-2t\right]\Biggr\rvert_1^2 = \frac{c}{6} \qquad \Longrightarrow \qquad c=6 \]

mypdf <- 

f <- function(t){
    ifelse((t>1 & t<2), (t-1)*(2-t), 0)
}

integral_c1 <- integrate(f, lower=1, upper=2)

cat('The constant c that guarantees the normalization is:', 1/integral_c1$value)

## The constant c that guarantees the normalization is: 6

b) Write the set of four R functions and plot the pdf and cdf, respectively

# Generic CDF
gen_cdf <- function(t, pdf, reltol = 1e-12){
  integrals = numeric(length(t))
  for (i in seq_along(t)){
    integrals[i] <- integrate(pdf, lower=-Inf, upper = t[i], rel.tol = reltol)$value
  }
  return(integrals)
}

# Generic inverse functions;
# bounds set such that it searches for a solution in [1, 2], the only interval in which the cdf is invertible
inverse <- function(f, lower=1, upper=2){
  function(y){
    uniroot(function(x){f(x) - y}, lower = lower, upper = upper)$root
  }
}

# Generic quantile function
gen_quantile <- function(y, cdf, lower.value=1, upper.value=2){
  output <- numeric(length(y))
  for (i in seq_along(y)){
    if (y[i]==1) output[i]<-upper.value       # boundary value
    else if (y[i]==0) output[i]<-lower.value  # boundary value
    else output[i] <- inverse(cdf)(y[i])
  }
  return(output)
}

# PDF
d_custom_1 <- function(t, c=6){
    ifelse((t>1 & t<2), c*(t-1)*(2-t), 0)
}

# CDF
p_custom_1 <- function(t){gen_cdf(t, d_custom_1)}

# Quantile function
q_custom_1 <- function(y){gen_quantile(y, p_custom_1)}

# Generate random numbers from the distribution
r_custom_1 <- function(n, seed=1){
  set.seed(seed)
  q_custom_1(runif(n))
}

Plot of the pdf and the cdf

time <- seq(0.5,2.5, by=0.05)

gg_distribution <- function(Y, Title='', xlabel='x', ylabel='y', col=TRUE){
  ggplot()+
    {
    if (col==TRUE) 
      geom_col(aes(x=factor(time), y=Y),  fill='gold', colour="black") 
    else
      geom_area(aes(x=factor(time), y=Y, group=1),  colour="black", fill = "lightblue")
    }+
  labs(title=Title, 
       x=xlabel,
       y=ylabel)+
  scale_x_discrete(breaks=seq(0.5, 2.5, by=0.1))+
  theme_bw()
}

gg_distribution(d_custom_1(time), "Probability density function", "Time (hours)", "f(x)")

gg_distribution(p_custom_1(time), "Cumulative density function", "Time (hours)", "F(x)")

I repeat the plot of the pdf and the cdf using the `highcharter’ package.

distributions_1 <- data.frame(X=factor(time, levels = time), pdf=d_custom_1(time), cdf=p_custom_1(time))

hc_pdf <- distributions_1 %>% hchart(
  'line', hcaes(x = time, y = pdf),
  color = "steelblue"
  ) %>%
  hc_title(text='Probability density function') %>% 
  hc_xAxis(title = list(text = 'Time (hours)')) %>%
  hc_yAxis(title = list(text = 'f(x)'))

hc_cdf <- distributions_1 %>% hchart(
  'line', hcaes(x = time, y = cdf),
  color = "steelblue"
  ) %>%
  hc_title(text='Cumulative density function') %>% 
  hc_xAxis(title = list(text = 'Time (hours)')) %>%
  hc_yAxis(title = list(text = 'F(x)'))

hc_pdf

hc_cdf

Test r_custom: sampling from a user’s defined pdf.

gg_sampling <- function(r_distribution, dp_distribution, nsamples, returnvalues = FALSE, Title='Sampling results from the defined pdf', xlabel='x', ylabel='Count (normalized)'){
  sampled = r_distribution(nsamples)
  gg <- ggplot() + 
    geom_histogram(aes(x=sampled, y=..density..),  binwidth = 0.05, center = 0.05, fill='gold', colour="black") +
    stat_function(fun=dp_distribution, color='firebrick') +
    labs(title=Title, 
        x=xlabel,
        y=ylabel)+
    scale_x_continuous(breaks=seq(0.7, 2.3, by=0.1), limits=c(0.7, 2.3)) +
    theme_bw()
  
  plot(gg)
  if (returnvalues) return(sampled)
}

gg_sampling(r_custom_1, d_custom_1, 3000, xlabel='Time (hours)')

c) Evaluate the probability that the student will finish the aptitude test in more than 75 minutes. And that it will take 90 and 120 minutes.

time_1 <- 75./60
time_2 <- 90./60
time_3 <- 120./60

cat('Probability that the student will finish the test in more than 75 minutes:', (1-p_custom_1(time_1))*100, '%\n')

## Probability that the student will finish the test in more than 75 minutes: 84.375 %

cat('Probability that the student will finish the test between 90 and 120 minutes:', (1-p_custom_1(90./60))*100, '%')

## Probability that the student will finish the test between 90 and 120 minutes: 50 %

The probability that it will take exactly 90 and 120 minutes is zero, because the distribution is continuous and \(P[X=a]=0 \; \forall a\) in a continuous distribution.

Exercise 6

The lifetime of tires sold by an used tires shop is \(10^4\cdot x\) km, where \(x\) is a random variable following the distribution function \[ \begin{equation} f(t) = \begin{cases} 2/x^2 & 1<x<2\\ 0 & \text{otherwise} \end{cases} \end{equation} \]

a) Write the set of four R functions and plot the pdf and cdf, respectively

# PDF
d_custom_2 <- function(x){
    ifelse((x>1 & x<2), 2/x^2, 0)
}

# CDF
p_custom_2 <- function(x){gen_cdf(x, d_custom_2, reltol=1e-7)}

# Quantile function
q_custom_2 <- function(y){gen_quantile(y, p_custom_2)}

# Generate random numbers from the distribution
r_custom_2 <- function(n, seed=1){
  set.seed(seed)
  q_custom_2(runif(n))
}

gg_distribution(d_custom_2(time), "Probability density function", "x (10^4 km)", "f(x)")

gg_distribution(p_custom_2(time), "Cumulative density function", "x (10^4 km)", "F(x)")

I repeat the plot of the pdf and the cdf using the `highcharter’ package.

distributions_2 <- data.frame(X=factor(time, levels = time), pdf=d_custom_2(time), cdf=p_custom_2(time))

hc_pdf <- distributions_2 %>% hchart(
  'line', hcaes(x = time, y = pdf),
  color = "steelblue"
  ) %>%
  hc_title(text='Probability density function') %>% 
  hc_xAxis(title = list(text = 'x')) %>%
  hc_yAxis(title = list(text = 'f(x)'))

hc_cdf <- distributions_2 %>% hchart(
  'line', hcaes(x = time, y = cdf),
  color = "steelblue"
  ) %>%
  hc_title(text='Cumulative density function') %>% 
  hc_xAxis(title = list(text = 'x')) %>%
  hc_yAxis(title = list(text = 'F(x)'))

hc_pdf

hc_cdf

b) Determine the probability that tires will last less than 15000 km

lifetime = 15000
x_test = lifetime/10^4
cat('The probability that tires will last less than 15000 km is', p_custom_2(x_test)*100, '%')

## The probability that tires will last less than 15000 km is 66.667 %

c) Sample 3000 random variables from the distribution and determine the mean value and the variance, using the expression \(Var(X) = E[X^2]-E[X]^2\)

# sample and plot at the same time using the function defined above
r_samples <- gg_sampling(r_custom_2, d_custom_2, 3000, returnvalues = TRUE, xlabel='x (10^4 km)')

# mean value and variance
exp_x <- mean(r_samples)
exp_x2 <- mean(r_samples^2)
var <- exp_x2-exp_x^2

cat('Mean value of x:', exp_x)

## Mean value of x: 1.3813

cat('\nVariance of x:', var)

## 
## Variance of x: 0.080035

glue('The average lifetime of the tires is {format(exp_x*10^4, digits=5)} km +- {format(sqrt(var)*10^4, digits=5)} km.')

## The average lifetime of the tires is 13813 km +- 2829 km.

Exercises Laboratory Session 02

Nicola Zomer

16 April, 2022

Packages and functions

Exercise 1

Exercise 2

A) Simulate the waiting time for 50 people at the doctor’s office and plot the relative histogram

B) What is the probability that a person will wait for less than 10 minutes?

C) Evaluate the average waiting time from the simulated data and compare it with the expected value

D) What is the probability for waiting more than one hour before being received?

Exercise 3

Task

Solution

Exercise 4

Task

Solution

Exercise 5

a) Using the integrate() R function, determine the constant \(c\) (and verify it analytically)

b) Write the set of four R functions and plot the pdf and cdf, respectively

c) Evaluate the probability that the student will finish the aptitude test in more than 75 minutes. And that it will take 90 and 120 minutes.

Exercise 6

a) Write the set of four R functions and plot the pdf and cdf, respectively

b) Determine the probability that tires will last less than 15000 km

c) Sample 3000 random variables from the distribution and determine the mean value and the variance, using the expression \(Var(X) = E[X^2]-E[X]^2\)