De cumulatieve verdelingsfunctie schetsen
import numpy as np
import matplotlib.pyplot as plt
Het verband zien tussen een continue dichtheidsfunctie $\varphi_X(x)$ en de overeenkomstige cumulatieve verdelingsfunctie $\Phi_X(x)$ is niet alleen nuttig qua inzicht. Op het examen zal je soms ook een schets van de ene grafiek moeten maken, gegeven de andere. Dit kan dus in twee richtingen werken, maar in deze blogpost focussen we ons op het schetsen van $\Phi_X(x)$ als $\varphi_X(x)$ gegeven is.
Ter herinnering: de cumulatieve verdelingsfunctie geeft de kans weer dat een toevalsvariabele $X$ een waarde aanneemt die kleiner is dan of gelijk is aan een bepaalde waarde $x$. Met andere woorden: $\Phi_X(x) = P(X \leq x)$.
Voorbeeld 1: rechthoekige dichtheidsfunctie¶
Laten we beginnen met een rechthoekige dichtheidsfunctie die de volgende vorm heeft:
$$ \varphi_X(x) = \begin{cases} 0.25 & \mid 1 \leq x \leq 5 \\ 0 & \mid \text{elders} \end{cases} $$
Anders gezegd: $X \sim \mathcal U(1, 5)$.
Deze dichtheidsfunctie is constant over het interval $[1, 5]$ met een hoogte van $0.25$. Dit betekent dat de totale oppervlakte onder de rechthoek (niet onverwacht) gelijk is aan $0.25 \cdot (5 - 1) = 1$.
De cumulatieve verdelingsfunctie $\Phi_X(x)$ is zoals steeds de oppervlakte onder de dichtheidsfunctie $\varphi_X(x)$ vanaf $-\infty$ tot $x$. In symbolen: $\int_{-\infty}^{x} \varphi_X(t)dt$. (We gebruiken hier de integratievariabele $t$ aangezien $x$ al gebruikt wordt om de bovengrens aan te geven.) Zo komen we tot volgende vorm voor $\Phi_X(x)$:
- Voor $x < 1$ is er geen kans, dus $\Phi_X(x) = 0$.
- Voor $1 \leq x \leq 5$ is de cumulatieve verdelingsfunctie gelijk aan de oppervlakte onder de constante $0.25$ van 1 tot $x$. De oppervlakte is dus gelijk aan $0.25 \cdot (x - 1)$.
- Voor $x > 5$ is de oppervlakte onder de curve gelijk aan 1, omdat de volledige oppervlakte onder de dichtheidsfunctie dan is bereikt.
De cumulatieve verdelingsfunctie $\Phi_X(x)$ is dus:
$$ \Phi_X(x) = \begin{cases} 0, & x < 1 \\ 0.25(x - 1), & 1 \leq x \leq 5 \\ 1, & x > 5 \end{cases} $$
def pdf_rect(x):
return 0.25 if 1 <= x <= 5 else 0
def cdf_rect(x):
if x < 1:
return 0
elif 1 <= x <= 5:
return 0.25 * (x - 1)
else:
return 1
x_vals = np.linspace(0, 6, 500)
pdf_vals = [pdf_rect(x) for x in x_vals]
cdf_vals = [cdf_rect(x) for x in x_vals]
plt.figure(figsize=(12, 10))
plt.subplot(2, 2, 1)
plt.title(r'Dichtheidsfunctie $\varphi_X(x)$')
plt.plot(x_vals, pdf_vals)
plt.fill_between(x_vals, pdf_vals, where=[1 <= x <= 5 for x in x_vals], alpha=0.1)
plt.xlabel('$x$')
plt.xlim(0, 6)
plt.ylim(0, 1.15)
plt.yticks(np.arange(0, 1.05, 0.25))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 2)
plt.title(r'Cumulatieve verdelingsfunctie $\Phi_X(x)$')
plt.plot(x_vals, cdf_vals)
plt.xlabel('$x$')
plt.xlim(0, 6)
plt.ylim(0, 1.15)
plt.yticks(np.arange(0, 1.05, 0.25))
plt.grid(linestyle='--', alpha=0.7)
bar_width = 0.5
bar_gap = 0.02
x_riemann = [1 + i * bar_width for i in range(8)]
pdf_riemann = [pdf_rect(x) for x in x_riemann]
cdf_riemann = np.cumsum([p * bar_width for p in pdf_riemann])
colors = [plt.cm.Blues(1 - i / len(x_riemann)) for i in range(len(x_riemann))]
plt.subplot(2, 2, 3)
plt.title(r'Gediscretiseerde dichtheidsfunctie $\pi_X(x)$')
# plt.plot(x_vals, pdf_vals)
# plt.fill_between(x_vals, pdf_vals, where=[1 <= x <= 5 for x in x_vals], alpha=0.1)
for i, x in enumerate(x_riemann):
plt.bar(x, pdf_riemann[i], width=bar_width-bar_gap, align='edge', color=colors[i])
plt.xlabel('$x$')
plt.xlim(0, 6)
plt.ylim(0, 1.15)
plt.yticks(np.arange(0, 1.05, 0.25))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 4)
plt.title(r'Gediscretiseerde cumulatieve verdelingsfunctie $\Phi_X(x)$')
# plt.plot(x_vals, cdf_vals)
for i in range(len(x_riemann)):
for j in range(i + 1):
plt.bar(
x_riemann[i],
pdf_riemann[j] * bar_width,
width=bar_width - bar_gap,
align='edge',
color=colors[j],
bottom=np.sum([p * bar_width for p in pdf_riemann[:j]]) if j > 0 else 0,
)
for x_repeat in [5, 5.5]:
for i in range(len(x_riemann)):
plt.bar(
x_repeat,
pdf_riemann[i] * bar_width,
width=bar_width - bar_gap,
align='edge',
color=colors[i],
bottom=np.sum([p * bar_width for p in pdf_riemann[:i]]) if i > 0 else 0,
)
plt.xlabel('$x$')
plt.xlim(0, 6)
plt.ylim(0, 1.15)
plt.yticks(np.arange(0, 1.05, 0.25))
plt.grid(linestyle='--', alpha=0.7)
Conclusie: een horizontaal lijnstuk in $\varphi_X(x)$ wordt een stijgende rechte in $\Phi_X(x)$. De rekenregel die hierachter zit is $\int h dx = hx + C$. Hoe hoger de horizontale lijn $y = h$ ligt, hoe sneller de rechte $hx$ dus zal stijgen. Merk ook op dat de rechte nooit kan dalen. Daarvoor zouden we een $h < 0$ nodig hebben, maar $\varphi_X(x) \geq 0$ dus dat kan in deze context niet.
Wie liever niet opnieuw nadenkt over de rekenregels rond integralen kan op de twee onderste grafieken ook intuitief zien waarom een horizontaal lijnstuk plots een stijgende rechte wordt. Als we de dichtheidsfunctie benaderen met een discrete kansmassafunctie $\pi_X(x) = 0.25$ voor $x \in \{1, 1.5, \ldots, 4.5\}$ en we berekenen daar de cumulatieve variant van, zien we hetzelfde patroon ontstaan.
Nu we dit weten, moeten we in de toekomst voor rechthoeken gewoon de totale oppervlakte bij het begin en bij het einde berekenen. In ons geval is de oppervlakte $P(X \leq 1) = 0.0$ en $P(X \leq 5) = 1.0$. We kunnen dus om $\Phi_X(x)$ te schetsen de twee coordinaten $(1, 0.0)$ en $(5, 1.0)$ verbinden met een rechte en we zijn klaar.
Helaas gaat het niet altijd zo eenvoudig zijn. Driehoekige functies bijvoorbeeld zullen geen rechte meer opleveren.
Voorbeeld 2: stijgende driehoekige dichtheidsfunctie¶
In dit voorbeeld hebben we een dichtheidsfunctie die de vorm heeft van een driehoek. Stel dat de dichtheidsfunctie wordt gegeven door:
$$ \varphi_X(x) = \begin{cases} 2x & \mid 0 \leq x \leq 1 \\ 0 & \mid \text{elders} \end{cases} $$
De cumulatieve verdelingsfunctie $\Phi_X(x)$ is opnieuw de oppervlakte onder de curve van de dichtheidsfunctie vanaf $-\infty$ tot $x$. In dit geval stijgt de oppervlakte niet lineair maar kwadratisch aangezien $\int (ax+b) dx = a \frac{x^2}{2} + bx + C$. Specifiek ik ons geval is $a=2$ en $b=0$, dus krijgen we $x^2 + C$. De cumulatieve verdelingsfunctie wordt dus (voor $C=0$):
$$ \Phi_X(x) = \begin{cases} 0 & \mid x < 0 \\ x^2 & \mid 0 \leq x \leq 1 \\ 1 & \mid x > 1 \end{cases} $$
De grafieken van $\varphi_X(x)$ en $\Phi_X(x)$ vind je hieronder. Het gediscretiseerde voorbeeld eronder toont opnieuw visueel aan hoe we van een stijgende rechte naar een parabool gaan zonder rekenregels voor integralen.
def pdf_triangle(x):
return 2 * x if 0 <= x <= 1 else 0
def cdf_triangle(x):
if x < 0:
return 0
elif 0 <= x <= 1:
return x**2
else:
return 1
x_vals = np.linspace(-0.1, 1.1, 500)
pdf_vals = [pdf_triangle(x) for x in x_vals]
cdf_vals = [cdf_triangle(x) for x in x_vals]
plt.figure(figsize=(12, 10))
plt.subplot(2, 2, 1)
plt.title(r'Dichtheidsfunctie $\varphi_X(x)$')
plt.plot(x_vals, pdf_vals)
plt.fill_between(x_vals, pdf_vals, where=[0 <= x <= 1 for x in x_vals], alpha=0.1)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 2)
plt.title(r'Cumulatieve verdelingsfunctie $\Phi_X(x)$')
plt.plot(x_vals, cdf_vals)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
bar_width_triangle = 0.1
bar_gap = 0.005
x_riemann_triangle = [i * bar_width_triangle for i in range(10)]
pdf_riemann_triangle = [pdf_triangle(x+0.05) for x in x_riemann_triangle]
cdf_riemann_triangle = np.cumsum([p * bar_width_triangle for p in pdf_riemann_triangle])
colors_triangle = [plt.cm.Blues(1 - i / len(x_riemann_triangle)) for i in range(len(x_riemann_triangle))]
plt.subplot(2, 2, 3)
plt.title(r'Gediscretiseerde dichtheidsfunctie $\pi_X(x)$')
# plt.plot(x_vals, pdf_vals)
# plt.fill_between(x_vals, pdf_vals, where=[0 <= x <= 1 for x in x_vals], alpha=0.1)
for i, x in enumerate(x_riemann_triangle):
plt.bar(x, pdf_riemann_triangle[i], width=bar_width_triangle-bar_gap, align='edge', color=colors_triangle[i])
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 4)
plt.title(r'Gediscretiseerde cumulatieve verdelingsfunctie $\Phi_X(x)$')
# plt.plot(x_vals, cdf_vals)
for i in range(len(x_riemann_triangle)):
for j in range(i + 1):
plt.bar(
x_riemann_triangle[i],
pdf_riemann_triangle[j] * bar_width_triangle,
width=bar_width_triangle - bar_gap,
align='edge',
color=colors_triangle[j],
bottom=np.sum([p * bar_width_triangle for p in pdf_riemann_triangle[:j]]) if j > 0 else 0,
)
for x_repeat in [1]:
for i in range(len(x_riemann_triangle)):
plt.bar(
x_repeat,
pdf_riemann_triangle[i] * bar_width_triangle,
width=bar_width_triangle - bar_gap,
align='edge',
color=colors_triangle[i],
bottom=np.sum([p * bar_width_triangle for p in pdf_riemann_triangle[:i]]) if i > 0 else 0,
)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
Net zoals voorheen kunnen we dus starten met het aanduiden van de twee coordinaten: in dit geval $(0, 0.0)$ en $(1, 1.0)$. Maar in plaats van de twee punten te verbinden met een rechte, moeten we ze nu verbinden met een dalparabool ($ax^2+bx+c$ met $a>0$).
Als de driehoek omgekeerd had gestaan (dus gevormd door een dalende rechte), dan moesten we de twee punten verbinden door een bergparabool ($ax^2+bx+c$ met $a<0$).
Voorbeeld 3: dalende driehoekige dichtheidsfunctie¶
In dit voorbeeld hebben we een dichtheidsfunctie die de vorm heeft van een driehoek die afneemt van $2$ bij $x = 0$ naar $0$ bij $x = 1$. De dichtheidsfunctie wordt als volgt gegeven:
$$ \varphi_X(x) = \begin{cases} 2(1 - x) & \mid 0 \leq x \leq 1 \\ 0 & \mid \text{elders} \end{cases} $$
Dit geeft ons de cumulatieve verdelingsfunctie:
$$ \Phi_X(x) = \begin{cases} 0 & \mid x < 0 \\ 2x - x^2 & \mid 0 \leq x \leq 1 \\ 1 & \mid x > 1 \end{cases} $$
Deze cumulatieve verdelingsfunctie heeft inderdaad de vorm van een bergparabool aangezien $a = -1 < 0$. Verder is de werkwijze helemaal identiek aan die van voorbeeld 2.
def pdf_reverse_triangle(x):
return 2 * (1 - x) if 0 <= x <= 1 else 0
def cdf_reverse_triangle(x):
if x < 0:
return 0
elif 0 <= x <= 1:
return 2 * x - x**2
else:
return 1
x_vals = np.linspace(-0.1, 1.1, 500)
pdf_vals = [pdf_reverse_triangle(x) for x in x_vals]
cdf_vals = [cdf_reverse_triangle(x) for x in x_vals]
plt.figure(figsize=(12, 10))
plt.subplot(2, 2, 1)
plt.title(r'Dichtheidsfunctie $\varphi_X(x)$')
plt.plot(x_vals, pdf_vals)
plt.fill_between(x_vals, pdf_vals, where=[0 <= x <= 1 for x in x_vals], alpha=0.1)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 2)
plt.title(r'Cumulatieve verdelingsfunctie $\Phi_X(x)$')
plt.plot(x_vals, cdf_vals)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
bar_width_triangle = 0.1
bar_gap = 0.005
x_riemann_triangle = [i * bar_width_triangle for i in range(10)]
pdf_riemann_triangle = [pdf_reverse_triangle(x + 0.05) for x in x_riemann_triangle]
cdf_riemann_triangle = np.cumsum([p * bar_width_triangle for p in pdf_riemann_triangle])
colors_triangle = [plt.cm.Blues(1 - i / len(x_riemann_triangle)) for i in range(len(x_riemann_triangle))]
plt.subplot(2, 2, 3)
plt.title(r'Gediscretiseerde dichtheidsfunctie $\pi_X(x)$')
for i, x in enumerate(x_riemann_triangle):
plt.bar(x, pdf_riemann_triangle[i], width=bar_width_triangle-bar_gap, align='edge', color=colors_triangle[i])
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(2, 2, 4)
plt.title(r'Gediscretiseerde cumulatieve verdelingsfunctie $\Phi_X(x)$')
for i in range(len(x_riemann_triangle)):
for j in range(i + 1):
plt.bar(
x_riemann_triangle[i],
pdf_riemann_triangle[j] * bar_width_triangle,
width=bar_width_triangle - bar_gap,
align='edge',
color=colors_triangle[j],
bottom=np.sum([p * bar_width_triangle for p in pdf_riemann_triangle[:j]]) if j > 0 else 0,
)
for x_repeat in [1]:
for i in range(len(x_riemann_triangle)):
plt.bar(
x_repeat,
pdf_riemann_triangle[i] * bar_width_triangle,
width=bar_width_triangle - bar_gap,
align='edge',
color=colors_triangle[i],
bottom=np.sum([p * bar_width_triangle for p in pdf_riemann_triangle[:i]]) if i > 0 else 0,
)
plt.xlabel('$x$')
plt.xlim(-0.1, 1.1)
plt.ylim(0, 2.15)
plt.yticks(np.arange(0, 2.15, 0.2))
plt.grid(linestyle='--', alpha=0.7)
Voorbeeld 4: complexe dichtheidsfunctie¶
We hebben een samengestelde dichtheidsfunctie die bestaat uit een trapezium en een rechthoek:
- Voor $1 \leq x \leq 2$: Een stijgende lijn van $0$ naar $0.20$
- Dit segment moeten we vertalen in een dalparabool van $(1, 0)$ naar $(2, 0.10)$
- Voor $2 \leq x \leq 5$: Een horizontale lijn t.h.v. $0.20$
- Dit segment moeten we vertalen in een stijgende rechte van $(2, 0.10)$ naar $(5, 0.70)$
- Voor $5 \leq x \leq 6$: Een dalende lijn van $0.20$ naar $0$
- Dit segment moeten we vertalen in een bergparabool van $(5, 0.70)$ naar $(6, 0.80)$
- Voor $7 \leq x \leq 8$: Een horizontale lijn t.h.v. $0.20$
- Dit segment moeten we vertalen in een stijgende rechte van $(7, 0.80)$ naar $(8, 1.00)$
- De tussenliggende waarden zijn nul
- Deze segmenten moeten we vertalen in horizontale lijnen
We kunnen dit als volgt beschrijven:
$$ \varphi_X(x) = \begin{cases} 0.20(x - 1) & \mid 1 \leq x \leq 2 \\ 0.20 & \mid 2 < x \leq 5 \\ -0.20(x - 6) & \mid 5 < x \leq 6 \\ 0.20 & \mid 7 \leq x \leq 8 \\ 0 & \text{elders} \end{cases} $$
Voor de liefhebbers (dit moet je niet kunnen uitschrijven op het examen):
$$ \Phi_X(x) = \begin{cases} 0 & \mid x < 1 \\ 0.10(x - 1)^2 & \mid 1 \leq x \leq 2 \\ 0.10 + 0.20(x - 2) & \mid 2 < x \leq 5 \\ 0.70 + 0.10(1 - (6 - x)^2) & \mid 5 < x \leq 6 \\ 0.80 & \mid 6 < x \leq 7 \\ 0.80 + 0.20(x - 7) & \mid 7 < x \leq 8 \\ 1 & \mid x > 8 \end{cases} $$
De constante $C$ die steeds uit de berekening van een integraal rolt, kiezen we hier zodanig dat elk segment van de grafiek mooi aansluit op het vorige.
def pdf_composite(x):
if 1 <= x <= 2:
return (x - 1) / 5
elif 2 < x <= 5:
return 0.20
elif 5 < x <= 6:
return -(x - 6) / 5
elif 7 <= x <= 8:
return 0.20
else:
return 0
def cdf_composite(x):
if x < 1:
return 0
elif 1 <= x <= 2:
return (x - 1)**2 / 10
elif 2 < x <= 5:
return 0.1 + 0.20 * (x - 2)
elif 5 < x <= 6:
return 0.7 + (1 - (6 - x)**2) / 10
elif 6 < x <= 7:
return 0.8
elif 7 < x <= 8:
return 0.8 + 0.2 * (x - 7)
else:
return 1
pdf_segments = {
"0 <= x < 1": {"x": np.linspace(0, 1, 100), "color": "gray"},
"1 <= x < 2": {"x": np.linspace(1, 2, 100), "color": "blue"},
"2 <= x < 5": {"x": np.linspace(2, 5, 100), "color": "green"},
"5 <= x < 6": {"x": np.linspace(5, 6, 100), "color": "red"},
"6 <= x < 7": {"x": np.linspace(6, 6.99, 100), "color": "orange"},
"7 <= x < 8": {"x": np.linspace(7, 8, 100), "color": "purple"},
"8 <= x < 9": {"x": np.linspace(8.01, 9, 100), "color": "gray"},
}
cdf_segments = {
"0 <= x < 1": {"x": np.linspace(0, 1, 100), "color": "gray"},
"1 <= x < 2": {"x": np.linspace(1, 2, 100), "color": "blue"},
"2 <= x < 5": {"x": np.linspace(2, 5, 100), "color": "green"},
"5 <= x < 6": {"x": np.linspace(5, 6, 100), "color": "red"},
"6 <= x < 7": {"x": np.linspace(6, 7, 100), "color": "orange"},
"7 <= x < 8": {"x": np.linspace(7, 8, 100), "color": "purple"},
"8 <= x < 9": {"x": np.linspace(8, 9, 100), "color": "gray"},
}
colors = [plt.cm.Blues(i / len(pdf_segments)) for i in range(len(pdf_segments))]
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.title(r'Dichtheidsfunctie $\varphi_X(x)$')
for label, segment in pdf_segments.items():
x_segment = segment["x"]
pdf_segment_vals = [pdf_composite(x) for x in x_segment]
plt.plot(x_segment, pdf_segment_vals, color=segment["color"])
plt.fill_between(x_segment, pdf_segment_vals, alpha=0.1, color=segment["color"])
plt.xlabel('$x$')
plt.xlim(0, 9)
plt.ylim(0, 1.05)
plt.yticks(np.arange(0, 1.05, 0.1))
plt.grid(linestyle='--', alpha=0.7)
plt.subplot(1, 2, 2)
plt.title(r'Cumulatieve verdelingsfunctie $\Phi_X(x)$')
for label, segment in cdf_segments.items():
x_segment = segment["x"]
cdf_segment_vals = [cdf_composite(x) for x in x_segment]
plt.plot(x_segment, cdf_segment_vals, color=segment["color"])
plt.xlabel('$x$')
plt.xlim(0, 9)
plt.ylim(0, 1.05)
plt.yticks(np.arange(0, 1.05, 0.1))
plt.grid(linestyle='--', alpha=0.7)
Conclusie¶
Om een schets te maken van $\Phi_X(x)$ kan je $\varphi_X(x)$ segment per segment vertalen. Eerst bepaal je de begin en eindcoordinaat o.b.v. de oppervlakte onder $\varphi_X(x)$ tot respectievelijk het begin en het einde van het segment. Vervolgens bepaal je de vorm van de lijn tussen de twee coordinaten die je nodig hebt:
| $\varphi_X(x)$ segment | $\Phi_X(x)$ segment |
|---|---|
| $0$ | horizontale lijn |
| horizontale lijn | stijgende rechte |
| stijgende rechte | dalparabool |
| dalende rechte | bergparabool |
Dubbelcheck op het einde steeds dat je op oppervlakte $1.0$ uitkomt.
