knitr::opts_chunk$set(dev = "ragg_png")
options(digits=5) # set number of digits equal to 5
# 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')
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)
The waiting time, in minutes, at the doctor’s is about 30 minutes, and the distribution follows an exponential pdf with rate 1/30.
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()
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 %
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
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 %
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.
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 %
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.
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.
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 %
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.
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
# 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)')
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.
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} \]
# 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
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 %
# 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.