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ě tři funkce:

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í

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.1 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.1.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
    1. Nahrajte data.  
    2. Vytvořte vektor pořadí hodnot.  
    3. Spočítejte pravděpodobnosti.  
    4. Vyneste do grafu.  
    5. Proveďte s pomocí funkce ecdf().

7.1.2 Kvantilová funkce

Kvantilová funkce je inverzní funkcí k 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.841621

7.1.3 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 vygenerovanou sekvencí.

Kód 
seed <- .Random.seed
head(seed, 10)
 [1]       10403           2  1194467966 -1032711624   835003675  -380439943
 [7]  -244224312  1955718902 -1402540239 -1417164957
Kó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.6957853
[1] 0.6957853 0.9502328

7.2 Diskrétní rozdělení

7.2.1 Poissonovo rozdělení

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
  1. Vytvořte jednoduchý simulátor deště pro jeden rok podle následujícího zadání
    a) Buď prší, nebo ne. Prvděpodobnost že neprší, je 90 %.
    b) Pokud pršelo předchozí časový úsek, tak se zvyšuje pravděpodobnost že prší dnes na 25 %.

Ř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

Pravděpodobnost, že \(x\) je menší nebo rovno 2.

Pravděpodobnost, že \(x\) je menší nebo rovno 2.

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.9540539

Hledaná 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
    1. 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}} \]
    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 > 12,4)=\)  
    3. Pro jakou hodnotu \(\gamma\) je pravděpodobnost přibližně:
      • \(P(X < \gamma) = 0,76\quad\gamma=\)  
      • \(P(X > \gamma) = 0,16\quad\gamma=\)

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 9.

Kód 
qt(p = 0.95, df = 100)
[1] 1.660234

7.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

ČHMÚ používá pro odhad a konstrukci \(m-\)denních průtoků 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 9) 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í
  1. Generujte 10 čísel s pomocí normálního rozdělení s parametry \(\mu = 2.3\) a \(\sigma = 4.23\).
  2. Generujte 50 čísel \(X\sim\mathsf{Po}(\lambda = 2.3)\) s pomocí set.seed(123). Doplňte střední hodnotu generovaného souboru  
  3. Jaká je pravděpodobnost, že veličina \(X\sim\mathsf{N}(1.3, 4)\), bude nabývat hodnot menších než 6?
  4. 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 ))
  5. Najděte následující hodnoty: - t_{0,95}(10), {0,1}^2(8), F{0,01}(13;4)