Probability and Statistics EHTP - HWRE AND CE

École Hassania des Travaux Publics

Department of Probability and statistics

Assignment Report
Practical work sheet
Students:
Otmane El Manssouri
Abdelkarim Fadel
Souhail El-Bouzaidi Drissi Date:
December 2025
Lecturer:
Mr. Hanini Mohammed
[Link]@[Link]
1
Probability and Statistics EHTP - HWRE AND CE

Exercise 1: Justification of the Inversion Method

for Simulation
1.1 Problem Statement
Let X be a discrete random variable taking values in
X(Ω) = {x1 , x2 , . . . , xn },
with probability mass function
n
X
P (X = xi ) = pi , pi ∈ (0, 1], pi = 1.
i=1

We define the cumulative sums:

i
X
si = pj , s0 = 0.
j=1

Since the sum of all probabilities is equal to 1, we have:

n
X
sn = pi = 1.
i=1

Let U be a random variable uniformly distributed on [0, 1], that is:

U ∼ U ([0, 1]).

We define the function:

n
X
g(u) = xi · 1(si−1 , si ] (u),
i=1

where 1A denotes the indicator function of the set A.

The simulated random variable is then defined as:
Y = g(U ).

The objective is to prove that Y has the same distribution as X.

2
Probability and Statistics EHTP - HWRE AND CE

1.2 Proof
The interval [0, 1] is divided into n subintervals determined by the cumulative
sums:
]s0 , s1 ], ]s1 , s2 ], . . . , ]sn−1 , sn ],
where the length of each subinterval is exactly pi , since:
X i  X i−1 
si − si−1 = pj − pj = pi .
j=1 j=1

The function g assigns:



x1 to any u ∈]s0 , s1 ],

x
2 to any u ∈]s1 , s2 ],
.
..


xn to any u ∈]sn−1 , sn ].


Thus, the value returned by g(U ) depends on the subinterval in which the
uniform random variable U falls.
We now compute P (Y = xi ).
From the definition of Y , we have:
Y = xi ⇐⇒ U ∈]si−1 , si ].

Therefore:
P (Y = xi ) = P (U ∈]si−1 , si ]).
Since U is uniformly distributed over [0, 1], the probability of falling in an
interval is equal to its length:
P (U ∈]si−1 , si ]) = si − si−1 .

But:
si − si−1 = pi .
Hence:
P (Y = xi ) = pi = P (X = xi ).
This equality holds for every i = 1, . . . , n. Therefore, the simulated variable
Y has the same distribution as X.

3
Probability and Statistics EHTP - HWRE AND CE

1.3 Application Python

import random as rd

def simul_X(values, probs):

# Calcul des bornes cumulées s_i
s = [0]
for p in probs:
[Link](s[-1] + p)

# Tirage uniforme
u = [Link]()

# Trouver l'intervalle où tombe u

for i in range(1, len(s)):
if s[i-1] < u <= s[i]:
return values[i-1]

Exercise 2: Discrete Random Variable Simula-

tion (Dart Game)
We consider a discrete random variable X modeling the score obtained from
a dart throw. The probabilities are:
P (X = 1) = 12 , P (X = 2) = 18 , P (X = 3) = 83 .

2.1 Calculation of Expectation and Variance

2.1.1 Expectation E[X]

E[X] = 1 · 12 + 2 · 18 + 3 · 3
8 = 4+2+9
8 = 15
8 = 1.875.

2.1.2 Variance Var(X)

First, we calculate the second moment E[X 2 ]:

E[X 2 ] = 12 · 21 + 22 · 81 + 32 · 3
8 = 4
8 + 84 + 27
8 = 35
8 = 4.375.

4
Probability and Statistics EHTP - HWRE AND CE

The variance is:

 2
Var(X) = E[X 2 ] − (E[X])2 = 35
8 − 15
8 = 280−225
64 = 55
64 ≈ 0.859375.

The standard deviation is:

p
σ= 55/64 ≈ 0.927.

2.2 Simulation Function simul X()

2.2.1 Sample Simulation and Empirical Mean X̄n

The function to simulate a realization of X, a sample of size n, and the

empirical mean are:

import random as rd

def simul_X():
U = [Link]()
if U < 0.5:
return 1
elif U < 0.5 + 0.125: # 0.625
return 2
else:
return 3

def simul_sample(n):
X = []
for i in range(n):
[Link](simul_X())
return X

def empirical_mean(X):
return sum(X) / len(X)

# Example for n = 10
X10 = simul_sample(10)
print("Sample X_10:", X10)
print("Mean X_10:", empirical_mean(X10))

5
Probability and Statistics EHTP - HWRE AND CE

2.3 Simulation Function simul Xn()

def simul_Xn(n, m):

"""
Simulate m empirical means, each based on a sample of size n.
Returns a list of m means.
"""
means = []
for k in range(m):
sample = simul_sample(n)
[Link](empirical_mean(sample))
return means

# Example usage:
means_10 = simul_Xn(10, 1)
print("Mean for n=10:", means_10[0])

2.3.1 Simulation Function simul Xn(10)

# Explicit simulation for n = 10 as in the report

X10 = [3, 1, 3, 1, 2, 2, 1, 3, 1, 1]
mean_10 = sum(X10) / len(X10)
print("X_10:", X10)
print("Empirical mean X_10:", mean_10)

2.4 Results and Convergence Interpretation

2.4.1 Calculation of the Sample Mean X̄10

The simulated sample for n = 10 is

x = [3, 1, 3, 1, 2, 2, 1, 3, 1, 1].

6
Probability and Statistics EHTP - HWRE AND CE

The empirical (sample) mean X̄10 is calculated as:

10
1 X 3 + 1 + 3 + 1 + 2 + 2 + 1 + 3 + 1 + 1 15
X̄10 = xi = = = 1.50.
10 i=1 10 10

The theoretical expectation is E[X] = 15/8 = 1.875. The distance |E[X] −

X̄10 | is then calculated:
|E[X] − X̄10 | = |1.875 − 1.50| = 0.375.

2.4.2 Simulated Results

The distance dn = |E[X] − X̄n | is calculated for increasing sample sizes n.

Size n Simulated Mean X̄n Distance |E[X] − X̄n |

10 1.50 0.375
100 1.76 0.115

2.5 Determination of Sample Size n (Question 5)

To determine the minimum sample size n such that the error |E[X] − X̄n | is
less than a given threshold ε > 0 with a maximum probability of error α, we
rely on the Bienaymé–Tchebychev Inequality. This inequality provides the
formal justification for the Law of Large Numbers.
The inequality states that the probability of the empirical mean deviating
from the expected value by more than ε is bounded by:

Var(X̄n )
P X̄n − E[X] ≥ ε ≤ .
ε2
Since the variance of the sample mean is Var(X̄n ) = Var(X)/n, we set the
maximum probability of error to α (e.g., α = 0.05):

Var(X)
P X̄n − E[X] ≥ ε ≤ ≤ α.
nε2
Solving for n, we obtain the required minimum sample size:
Var(X)
n≥ .
αε2
7
Probability and Statistics EHTP - HWRE AND CE

Numerical Example (α = 0.05, ε = 0.01): Using the calculated variance

Var(X) = 55/64 ≈ 0.859375:
0.859375 0.859375
n≥ 2
= ⇒ n ≥ 171875.
0.05 · (0.01) 0.000005

Therefore, a minimum sample size of n = 171,875 is required to guarantee

that the deviation is less than 0.01 with a probability of error no more than
5%.

Exercise 3
We consider the number X of screening tests performed on a given day, where
X ∼ Poisson(λ = 5).
Each test is independently positive with probability p = 0.1. Let Y denote
the number of positive tests.

1. Simulation of the variable Y

Conditionally on X = x, the number of positive tests is distributed as
Y | (X = x) ∼ Binomial(x, 0.1).

The following Python function simulates a realization of Y using this condi-

tional structure:

import [Link] as rd

def simulate_Y():
"""
Simulates the number of positive tests Y.
X ~ Poisson(5)
Y | X=x ~ Binomial(x, 0.1)
"""
X = [Link](5) # Number of tests performed
Y = [Link](X, 0.1) # Number of positive tests given X
return Y

8
Probability and Statistics EHTP - HWRE AND CE

2. Simulation of a sample and histogram comparison

Since Y theoretically follows a Poisson distribution with parameter p·λ = 0.5,
we compare two samples:

• one simulated using repeated calls of simulate Y(), and

• one drawn directly from a Poisson(0.5) distribution.

The following code generates the samples and produces the histogram com-
parison:

import numpy as np
import [Link] as plt
import [Link] as rd

n = 10000

# Sample using the simulation model

sample_simulated = [simulate_Y() for _ in range(n)]

# Sample directly from Poisson(0.5)

sample_poisson = [Link](0.5, size=n)

# Plot the histogram

[Link](figsize=(10,5))

[Link](sample_simulated, bins=15, density=True, alpha=0.6,

label="Simulated Y")

[Link](sample_poisson, bins=15, density=True, alpha=0.6,

label="Poisson(0.5)", color="red")

[Link]("Number of positive tests")

[Link]("Relative frequency")
[Link]()
[Link]("Comparison between empirical distribution of Y and Poisson(0.5)")
[Link](True)
[Link]()

9
Probability and Statistics EHTP - HWRE AND CE

Figure 1: Empirical comparison between simulated Y and a Poisson(0.5) law.

We observe that the empirical histogram issued from the compound simula-
tion almost overlaps perfectly with the histogram generated from a theoretical
Poisson(0.5) distribution. This confirms numerically the thinning property of
the Poisson distribution, i.e., that if X ∼ Poisson(5) and each event is kept
independently with probability p = 0.1, then the resulting count of retained
events follows Poisson(p · 5) = Poisson(0.5).

3. Justification
By the Poisson thinning theorem, if

X ∼ Poisson(λ) and Y | X = x ∼ Binomial(x, p),

then the marginal distribution of Y is

Y ∼ Poisson(pλ).

In our case,
p · λ = 0.1 × 5 = 0.5,
which matches the empirical distribution displayed in the histogram. Thus,
increasing the sample size confirms numerically that the simulation is consis-
tent with the theoretical model.

10
Probability and Statistics EHTP - HWRE AND CE

Exercise 4
A geometric distribution models the number of independent Bernoulli trials
required to obtain the first success. Let X denote this waiting time, where
success occurs with probability p at each trial. Therefore,

P (X = k) = (1 − p) k−1 p, for k = 1, 2, 3, . . .

1. Simulation of a geometric random variable

To simulate such a variable, we repeatedly generate independent Bernoulli
trials until the first success occurs. We use [Link](), which generates a
number uniformly in [0, 1]. If the generated number is less than p, we consider
a success; otherwise, the trial is repeated.
The following Python function implements this simulation:

import [Link] as rd

def Geometrique(p):
"""
Simulates a geometric random variable X ~ G(p)
representing the number of trials until the first success.

if not (0 < p <= 1):

raise ValueError("p must be in (0, 1]")

k = 1 # Starting at trial number 1

while True:
U = [Link]() # U ~ Uniform(0, 1)
if U < p: # Success condition
return k
k += 1 # Failure, continue

11
Probability and Statistics EHTP - HWRE AND CE

Exercise 5
1. Equivalence property of the pseudo-inverse function
Let F be a cumulative distribution function (CDF), and define its pseudo-
inverse (quantile function) by

F −1 (u) = inf{x ∈ R : F (x) ≥ u}, u ∈ [0, 1].

We show that
y ≤ F (x) ⇐⇒ F −1 (y) ≤ x.

Proof of (F −1 (y) ≤ x ⇒ y ≤ F (x))

Assume F −1 (y) ≤ x. By definition of the infimum, F −1 (y) is the smallest

value such that F (z) ≥ y. Therefore,

F (F −1 (y)) ≥ y.

Since F is non-decreasing, we have

F (x) ≥ F (F −1 (y)) ≥ y,

which implies
y ≤ F (x).

Proof of (y ≤ F (x) ⇒ F −1 (y) ≤ x)

Assume y ≤ F (x). Then x belongs to the set

A = {z ∈ R : F (z) ≥ y}.

Since F −1 (y) is the infimum of A, we must have

F −1 (y) ≤ x.

Combining the two implications, we conclude:

y ≤ F (x) ⇐⇒ F −1 (y) ≤ x.

12
Probability and Statistics EHTP - HWRE AND CE

Furthermore, equality holds when y is exactly in the range of F , i.e., when

there exists some x0 such that F (x0 ) = y. In this case, F −1 (y) is the smallest
such x0 .

2. Distribution of the variable X = F −1(U )

Let U ∼ Uniform(0, 1) and define
X = F −1 (U ).
Then we compute the CDF of X:
P (X ≤ x) = P (F −1 (U ) ≤ x).
Using the equivalence proven above, we obtain
P (F −1 (U ) ≤ x) = P (U ≤ F (x)).
Since U is uniform over [0, 1], its CDF is:
P (U ≤ t) = t, t ∈ [0, 1].
As F (x) ∈ [0, 1], it follows that
P (X ≤ x) = F (x).

Hence, the random variable X = F −1 (U ) has cumulative distribution function

F . This fundamental result is the basis of the inverse transform method used
in probability simulation.

Exercise 6
We aim to simulate samples from an exponential distribution using the in-
version method.

1. Inverse transform sampling for an exponential law

The cumulative distribution function (CDF) of an exponential random vari-
able with parameter λ > 0 is
F (x) = 1 − e−λx , x ≥ 0.

13
Probability and Statistics EHTP - HWRE AND CE

Let U ∼ Uniform(0, 1). Setting u = F (x), we obtain

1
u = 1 − e−λx ⇒ x = − ln(1 − u).
λ

Since 1 − u is also uniform on [0, 1], the implementation may equivalently use
1
X = − ln(U ).
λ
This constitutes the basis of the inverse transform method: starting from a
uniform sample, we generate exponential realizations.

2. Python implementation
The following Python function implements inverse transform sampling for
n independent exponential samples. This is the code extracted from our
implementation:

import numpy as np
import [Link] as plt

def inverse_transform_exponential(lambda_param, n_samples):

U = [Link](0, 1, n_samples)
X = -(1.0 / lambda_param) * [Link](U)
return X

# --- Simulation Parameters ---

LAMBDA = 0.5
N_SAMPLES = 10000

# --- Generate Samples ---

simulated_samples = inverse_transform_exponential(LAMBDA, N_SAMPLES)

# --- Visualization ---

# 1. Plot a Histogram of the Simulated Samples

[Link](figsize=(10, 6))
[Link](simulated_samples, bins=50, density=True, alpha=0.6,
,→ color='skyblue',

14
Probability and Statistics EHTP - HWRE AND CE

label='Simulated Samples')

# 2. Plot the Theoretical PDF for Comparison

x_range = [Link](0, [Link](simulated_samples), 500)
theoretical_pdf = LAMBDA * [Link](-LAMBDA * x_range)
[Link](x_range, theoretical_pdf, 'r-', linewidth=2,
label=f'Theoretical PDF ($\\lambda={ LAMBDA} $)')

[Link](f'Inverse Transform Sampling for Exponential Distribution

,→ (N={ N_SAMPLES} )')
[Link]('X (Simulated Value)')
[Link]('Density')
[Link]()
[Link](axis='y', linestyle='--')
[Link](left=0)
[Link]()

Figure 2: Inverse transform sampling for exponential distribution (N=1000).

3. Interpretation of results
The generated histogram closely follows the red theoretical exponential curve
derived from the PDF
f (x) = λe−λx .

15
Probability and Statistics EHTP - HWRE AND CE

This visual match confirms that:

• the inverse transform sampling has produced a correct exponential sam-

ple,

• the empirical density converges toward the theoretical density as the

number of samples increases.

When λ = 0.5, the theoretical mean is

1
E[X] = = 2.
λ

Numerically, calculating the empirical average of the simulated values gives

a value close to 2, which confirms the theoretical result and validates the
correctness of the implementation.

16