Kód
r___() # generování náhodných čísel z rozdělení
d___() # funkce hustoty rozdělení
p___() # Pravděpodobnostní funkce
q___() # Kvantilová funkce rozdělení7 Rozdělení pravděpodobnosti
Cíle cvičení
- Vědět kde v R najít rozdělení.
- Umět vysvětlit vztah mezi jednotlivými funkcemi rozdělení.
- Umět počítat s náhodnými čísly z rozdělení.
V software R jsou všechna základní rozdělení součástí balíku stats. A jejich výčet získáme vyvoláním nápovědy ?distributions. Pracovat s rozděleními budeme čtyřmi způsoby: - budeme generovat náhodná čísla z rozdělení - budeme pracovat s kvantilovou funkcí daného rozdělení - počítat s hustotou a pravděpododobnostní funkcí
Balík stats obsahuje mnoho rozdělení. Alternativní, Binomické, Multinomické, Poissonovo, Normální, Fischer-Snedecorovo, \(\chi^2\), Studentovo \(t\), Gumbellovo a mnohá další.
Podívejte se na ?Distributions.
7.1 Funkce rozdělení
Pro práci s rozděleními jsou k dispozici obecně čtyři funkce:
7.1.1 Hustota pravděpodobnosti
Probability density function (PDF)
Hustota pravděpodobnosti je funkce, které popisuje rozložení pravděpodobnosti u spojitě rozdělené náhodné veličiny, tz. vyjadřuje, jak je pravděpodobnost rozprostřena kolem dané hodnoty \(x\).
- Pravděpodobnost je nezáporná \[f(x)>0;\forall x\]
- Celková pravděpodobnost je rovna 1 \[\int\limits_{-\infty}^\infty f(x)dx=1\]
- Pravděpodobnost, že veličina spadne do intervalu \([a,b]\) je \[P(a \leq X \leq b) = \int\limits_a^bf(x)dx\]
- Pro spojité rozdělení platí \[P(X=x)=0\]
Pro úplnost jsou zde uvedeny funkce pro permutace factorial() a kombinace choose().
Kód
par(mfrow = c(1,2))
plot(x = 0:6,
y = factorial(0:6),
type = "s",
main = "Permutace",
ylab = bquote("log(0!:6!)"),
xlab = "",
log = "y")
plot(x = 0:6,
y = choose(6, 0:6),
type = "s",
xlab = "",
main = "Binomické koeficienty")- 1
- Kolika způsoby je možné seřadit \(n\leq6\) prvků?
- 2
- Kolika způsoby lze vybrat \(n-\)tici z \(6\) prvků?
7.1.2 Distribuční funkce rozdělení
Rozdělení pravděpodobnosti náhodné veličiny lze jednoznačně popsat tzv. distribuční funkcí Rovnice 7.1.
\[ F(x) = P(X\leq x) = P(\omega_i\in\Omega : X(\omega_i)\leq x) \tag{7.1}\]
7.1.2.1 Empirická distribuční funkce
Sestavuje se na základě sesbíraného vzorku dat, sděluje jak velký podíl dat je stejných nebo menších než konkrétní hodnota. K jejímu vytvoření potřebujeme data seřazená.
Kód
Qd <- read.fwf(file = "./data/01138000.dly",
widths = c(8, rep(10, times = 5)))[, 5]\[ P(X\leq x_k) \approx \dfrac{k - 0.3}{n + 0.4} \] kde \(k\) je pořadí v seřazeném souboru a \(n\) je počet prvků v souboru.
Úloha
- Nahrajte data.
- Vytvořte vektor pořadí hodnot.
- Spočítejte pravděpodobnosti.
- Vyneste do grafu.
- Proveďte s pomocí funkce
ecdf().
- Nahrajte data.
7.1.3 Kvantilová funkce
Kvantilová funkce je inverzní k funkci distribuční. Pokud distribuční funkce udává s jakou pravděpodobností bude hodnota náhodného pokusu menší nebo rovna \(x\), kvantilová funkce udává, pro jaké \(x\) bude výsledek náhodného pokusu s danou pravděpodobností \(y\) měnší nebo roven \(x\).
Kód
qnorm(p = 0.8, mean = 2, sd = 1)- 1
- S \(80\%\) pravděpodobností bude hodnota $$2.8416212.
[1] 2.8416217.1.4 Náhodná čísla z rozdělení
Pro generování náhodných čísel lze použít rozdělení.
Kód
runif(n = 10, min = 0, max = 1)
rpois(n = 15, lambda = 2.4)- 1
- Generování 10 čísel z rovnoměrného rozdělení z intervalu \((0;1)\)
- 2
- Generování 15 čísel z Poisonova rozdělení z intervalu \((0;1)\)
Generovaná čísla nejsou náhodná v pravém slova smyslu, ale označují se jako pseudonáhodná, neboť při jejich tvorbě se vychází z jiste sekvence čísel. Tuto sekvenci je možné přímo zvolit, čímž je zajištěna funkcí set.seed() volanou před každou vytvořenou množinou čísel.
Kód
seed <- .Random.seed
head(seed, 10) [1] 10403 2 530665372 1642361759 -633923083 -1094247284
[7] -1229413022 -1745929923 1979961783 -360458658Kód
plot(.Random.seed)
Kód
x <- sample(x = 1:1e3, size = 1)
set.seed(x)
runif(1)
set.seed(x)
runif(2)[1] 0.9147465
[1] 0.9147465 0.8205029
Poznámka
| Rozdělení | Typ / vlastnosti | Typické použití | Funkce v R (balík stats) | Příklad použití |
|---|---|---|---|---|
| Normální (Gaussovo) | Symetrické, definované střední hodnotou a směrodatnou odchylkou | Měřící chyby, přirozená variabilita, CLT | dnorm, pnorm, qnorm, rnorm |
rnorm(1000, mean=0, sd=1) |
| Lognormální | Kládně vymezené, silně pravostranné | Intenzita srážek, ekonomické veličiny, multiplikativní procesy | dlnorm, plnorm, qlnorm, rlnorm |
rlnorm(1000, meanlog=0, sdlog=1) |
| Exponenciální | Paměťově bezstarostné, kladná osa | Čekací doby, doby mezi událostmi (Poissonův proces) | dexp, pexp, qexp, rexp |
rexp(1000, rate=0.2) |
| Gamma | Kladné, flexibilní tvar (shape–rate) | Modelování srážek, nízkých průtoků, hydrologické extrémy | dgamma, pgamma, qgamma, rgamma |
rgamma(1000, shape=2, rate=0.5) |
| Chí-kvadrát (χ²) | Kladné, zvláštní případ gamma | Testy rozptylu, LRT testy | dchisq, pchisq, qchisq, rchisq |
rchisq(1000, df=10) |
| Studentovo t | Těžké ocasy, symetrické | Odhady střední hodnoty při malých vzorcích | dt, pt, qt, rt |
rt(1000, df=5) |
| F-rozdělení | Poměr dvou rozptylů | ANOVA, testy poměru rozptylů | df, pf, qf, rf |
rf(1000, df1=5, df2=20) |
| Rovnoměrné (Uniform) | Konstantní hustota | Náhodné vzorkování v intervalu, neinformativní předpoklady | dunif, punif, qunif, runif |
runif(1000, 0, 1) |
| Weibullovo | Kladné, flexibilní | Spolehlivost, doby selhání, délky suchých období | dweibull, pweibull, qweibull, rweibull |
rweibull(1000, shape=2, scale=5) |
| Beta | Omezené na (0, 1) | Poměry, pravděpodobnosti, Bayesovské prior distribuce | dbeta, pbeta, qbeta, rbeta |
rbeta(1000, 2, 5) |
| Cauchyho | Extrémně těžké ocasy, bez střední hodnoty a rozptylu | Robustní modelování, odlehlé hodnoty | dcauchy, pcauchy, qcauchy, rcauchy |
rcauchy(1000) |
7.2 Diskrétní rozdělení
Diskrétní rozdělení jsou určena výčtem hodnot spolu s pravděpodobností jejich nabytí.
Poznámka
| Rozdělení | Typ / vlastnosti | Typické použití | Funkce v R (balík stats) | Příklad použití |
|---|---|---|---|---|
| Poissonovo | Počty událostí v čase/prostoru | Počty srážkových událostí, příchody, počet výskytů | dpois, ppois, qpois, rpois |
rpois(1000, lambda=3) |
| Binomické | Počet úspěchů v pevně daném počtu pokusů | Pravděpodobnost výskytu jevu (např. déšť/bez deště) | dbinom, pbinom, qbinom, rbinom |
rbinom(1000, size=10, prob=0.4) |
| Negativně binomické | Počet neúspěchů do dosažení k úspěchů, nebo overdisperzní Poisson | Nadpoissonovská variabilita srážek, počet epizod | dnbinom, pnbinom, qnbinom, rnbinom |
rnbinom(1000, size=5, mu=2) |
| Geometrické | Počet neúspěchů do prvního úspěchu | První výskyt jevu, jednoduché čekací doby | dgeom, pgeom, qgeom, rgeom |
rgeom(1000, prob=0.3) |
| Hypergeometrické | Výběr bez vracení | Losování, ekologické vzorkování, bez-náhrady | dhyper, phyper, qhyper, rhyper |
rhyper(1000, m=25, n=75, k=10) |
| Multinomické | Více kategorií, opakované pokusy | Simulace kategoriální volby, více-třídní jevy | rmultinom |
rmultinom(10, size=20, prob=c(.2,.5,.3)) |
7.2.1 Binomické rozdělení
Popisuje počet \(m\) úspěchů v posloupnosti \(N\) nezávislých 0/1 pokusů.
\[ X\sim\mathsf{Bi}(n,p) \] ::: callout-tip ## Úloha
- V prosinci je pravděpodobnost sněhově srážkového dne \(p=0.15\). Spočítejte pravděpodobnost, že počet takových dnů bude 20.
Kód
pbinom(q = 20/31, size = 31, prob = 0.15)[1] 0.006486146:::
7.2.2 Poissonovo rozdělení
Poissonovo rozdělení (jinak také rozdělení vzácných jevů) je získáno aproximací binomického rozdělení v případech, kdy počet pozorování se limitně blíží k nekonečnu (\(N\rightarrow \infty\)). Rozdělení popisuje počet událostí, v daném časovém nebo prostorovém intervalu, pokud tyto nastávají nezávisle a střední doba (intenzita) výskytu je konstantní.
7.2.2.1 Hustota pravděpodobnosti
kde \(k\) je počet výskytů jevu a \(\lambda\) je střední počet výskytů v intervalu. \[ P(X = k) = \dfrac{\lambda^ke^{-k}}{k!} \]
Úloha
Průměrně se vyskytuje během zimy 5 výrazných sněhových epizod za měsíc. Jaká je pravděpodobnost, že v daném zimním měsící se vyskytnou přesne 3 epizody?
\[ \begin{align} P(X = 3) = \dfrac{5^3e^{-5}}{3!}\\ P(X = 3) = \dfrac{125}{6}e^{-5}\approx 20,83\cdot 0,006738 = 0,1404 \end{align} \]
Pravděpodobnost, že teplotní senzor během měření selže během jakéhokoliv dne je
Kód
# par(mfrow = c(2, 2))
lambda <- c(1, 4, 8)
plot(x = 1:20, ylim = c(0,0.4), type = "n")
for(l in lambda) {
lines(dpois(x = 1:20, lambda = l), lty = l)
points(dpois(x = 1:20, lambda = l), col = l)
}
Úloha
- Vytvořte jednoduchý simulátor deště pro jeden rok podle následujícího zadání
- Buď prší, nebo ne. Prvděpodobnost že neprší, je 90 %.
- Pokud pršelo předchozí časový úsek, tak se zvyšuje pravděpodobnost že prší dnes na 25 %.
- Buď prší, nebo ne. Prvděpodobnost že neprší, je 90 %.
Řešení:
Kód
x <- vector(mode = "numeric", length = 365)
for(i in seq_along(x)) {
if(i == 1) next
x[i] <- ifelse(
test = x[i - 1] > 0,
yes = rbinom(n = 1,
size = 1,
prob = 0.10) * rchisq(n = 1,
df = 1,
ncp = 0),
no = rbinom(n = 1,
size = 1,
prob = 0.25) * rchisq(n = 1,
df = 3,
ncp = 1))
}
barplot(x,
ylim = 1.5*c(max(x), 0),
col = "dodgerblue4",
border = "dodgerblue4")- 1
- Předem alokuejeme vektor hodnot.
- 2
- Je třeba ošetřit případ, kdy začínáme a neexistuje úsek \(x(t-1)\).
- 3
- Pokud pršelo v předchozím čase
7.3 Spojitá rozdělení
Během počítání s rozděleními se snažíme nalézt odpovědi na podobné otázky:
“Jaká je pravděpodobnost, že veličina \(X\) je menší než 7.5?”
7.3.1 Normální rozdělení
Chceme spočítat pravědpodobnost, že náhodná hodnota, vybraná z populace \(X\sim \mathsf{N}(\mu = 3; \sigma^2 = 2,2)\) bude menší než \(5,5\). Početně je úkon řešen standardizací a určením nového
Provedeme standardizaci na \(Z\sim\mathsf{N}(0; 1)\) \[ Z = \dfrac{X - \mu}{\sigma} \approx \dfrac{5,5 - 3}{1,48324} \approx 1,6855 \] \(P(Z < 1,69)\) Spočítáme s pomocí kvantilu Standardizovaného Normálního rozdělení. Tento úkon R dělá z nás, ale pro porovnání si uvedeme obě varianty:
Kód
cbind(
N = pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE),
Z = pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE)
) N Z
[1,] 0.9540539 0.9540539Hledaná pravděpodobnost je 95.45 %.
Poznámka
Zde je dobré upozornit na numerické zaokrouhlovací chyby, které během počítačového zpracování dat vznikají.
Kód
all.equal(
pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE),
pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE)
)[1] "Mean relative difference: 3.473763e-08"Kód
all.equal(
pnorm(q = 1.6855, mean = 0, sd = 1, lower.tail = TRUE),
pnorm(q = 5.5, mean = 3, sd = sqrt(2.2), lower.tail = TRUE), tolerance = 0.005
)[1] TRUE
Úloha
- S pomocí funkce
curve()vykreslete hustotu funkce normovaného normálního rozdělení danou předpisem: \[ f(x) = \dfrac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} \] - Máme náhodnou veličinu \(X\), jež sleduje Normální rozdělení s parametry \(\mu = 5\) a \(\sigma^2 = 4,13\), nebo-li \(X\sim \mathsf{N}(5; 4,13)\). Určete pravděpodobnost:
- \(P(X > 8)=\)
- \(P(X \leq 8)=\)
- \(P(X > 26,4)=\)
- \(P(X > 8)=\)
- Pro jakou hodnotu \(\gamma\) je pravděpodobnost přibližně:
- \(P(X < \gamma) = 0,76\quad\gamma=\)
- \(P(X > \gamma) = 0,16\quad\gamma=\)
- \(P(X < \gamma) = 0,76\quad\gamma=\)
- S pomocí funkce
- V Kapitola 5 jsme si představovali funkci
curve(). Pro praktické počítání s Normální rozdělením je užitečné si pamatovat pravděpodobnosti vybraných intervalů \(\mu\pm\{1,2,3\}\sigma\). S pomocí funkceintegrate()
\[ \int^{\mu+2\sigma}_{\mu-2\sigma}\dfrac{1}{\sigma\sqrt{2\pi}}\exp\left(-\dfrac{(x-\mu)^2}{2\pi^2}\right)dx \approx 0.95 \]
- V Kapitola 5 jsme si představovali funkci
7.3.2 Studentovo-\(t\) rozdělení
Používá se pro odhad střední hodnoty souboru pocházejícího z normálního rozdělení s neznámým parametrem \(\sigma\). Více v Kapitola 10.
Kód
qt(p = 0.95, df = 100)[1] 1.6602347.3.3 Logaritmicko-normální rozdělení
Nabývá pouze kladných hodnot. Můžeme jej nalézt například v rozdělení **průtoků, simulacích chemické konncentrace.
Poznámka
Pro odhad a konstrukci \(m-\)denních průtoků používá Český hydrometeorologický ústav logaritmicko-normální rodzdělení s pěti parametry.
tvar transformace se třemi, s pěti a ve standardizovaném tvaru pro průtoky
- \(\mu\): výběrový průmer normálního rozdělení
- \(\sigma\): směrodatná odchylka normálního rozdělení
- \(\gamma\): parametru posunu
- \(\alpha\): parametr měřítka
- \(\beta\): parametr tvaru (3. a 4. moment)
\[ f(x;\mu, \sigma, \gamma, \alpha, \beta) = \dfrac{\beta}{(x-\gamma)\sigma\sqrt{2\pi}}\exp\left(-\frac{\left(\ln{\left(\dfrac{x-\gamma}{\alpha}\right)}-\mu\right)^2}{2\sigma^2}\right),\quad x>\gamma \]
7.3.4 \(F\) rozdělení
Jinak také Fisherovo-Snedecorovo rozdělení je používáno pro sestrojení \(100(1-\alpha)\%\) intevalu spolehlivosti pro podíl rozptylů normálního rozdělení a je to modelové rozdělení testovací statistiky pro ověření shodnosti dvou rozptylů. V ?sec-anova nalezneme \(F-\)rozdělení u testování hypotézy o rovnosti středních hodnot u více než dvou výběrových souborů.
7.3.5 \(\chi^2\) rozdělení
Narozdíl od předchozích zde uvedených rozdělení není \(\chi^2\) rozdělení u pozorovaných veličin příliš časté. Nicméně je to významné rozdělení z hlediska testování statistických hypotéz (Kapitola 10) při porovnávání rozdělení vzájemně. Používá se při stanovení intervalů spolehlivosti výběrový rozptyl.
\[ \dfrac{(n-1)s_x^2}{\chi^2_{1-\alpha/2}} \leq \sigma_x^2\leq \dfrac{(n-1)s_x^2}{\chi^2_{\alpha/2}} \tag{7.2}\]
Cvičení
- Generujte 10 čísel s pomocí normálního rozdělení s parametry \(\mu = 2.3\) a \(\sigma = 4.23\).
- Generujte 50 čísel \(X\sim\mathsf{Po}(\lambda = 2.3)\) s pomocí
set.seed(123). Doplňte střední hodnotu generovaného souboru
- Jaká je pravděpodobnost, že veličina \(X\sim\mathsf{N}(1.3, 4)\), bude nabývat hodnot menších než 6?
- Předpokládejme, že \(s_x^2=8,6\) při \(15\) stupních volnosti. Doplňte následující výraz k určení 95% intervalu spolehlivosti pro \(\sigma_x^2\) (viz Rovnice 7.2).
c(15*/qchisq(,),15*/qchisq(,15))
- Najděte následující hodnoty: - t_{0,95}(10), ^2(8), F*{0,01}(13;4)