Data obsahují počet narozených dětí podle pohlaví v České republice od roku 1947 do 2018, jsou dostupná na webu https://www.mortality.org/. Vyzkoušíme na nich tři příklady - testování hypotézy, test dobré shody a lineární regresi.

Příklad 1

Na několika internetových portálech se můžeme dočíst, že poměr narození chlapců a dívek je 105:100. Z našich dat chceme ověřit, jestli to skutečně platí.

p_chlapec = 105 / (105 + 100) # převedeme poměr ze zadání na pravděpodobnost
cat("Pravděpodobnost narození chlapce:", p_chlapec, "\n")
## Pravděpodobnost narození chlapce: 0.5121951
# H_0: p = p_chlapec (= 0.51)
# H_1: p != p_chlapec
# Provedeme dvoustranný test s hladinou významnosti 5%.

df = read.csv("CZEbirth.txt")
data = df$Births

n = length(data)/2
chlapci = n*0 # vektor počtu narozených chlapců pro každý rok
soucty = n*0 # vektor počtu narozených dětí pro každý rok

for (i in 1:n)
{
  chlapci[i] = data[i*2-1] # index násobíme dvěma, protože v souboru jsou střídavě chlapci a dívky
  soucty[i] = data[i*2-1] + data[i*2]
}

pravdepodobnosti = chlapci / soucty

vyberovy_prumer = mean(pravdepodobnosti) # ekvivalentně vyberovy_prumer = sum(pravdepodobnosti) / n
cat("Výběrový průměr:", vyberovy_prumer, "\n")
## Výběrový průměr: 0.5139382
vyberovy_rozptyl = var(pravdepodobnosti) # ekvivalentně vyberovy_rozptyl = sum((pravdepodobnosti-vyberovy_prumer)^2)/(n-1)
cat("Výběrový rozptyl:", vyberovy_rozptyl, "\n\n")
## Výběrový rozptyl: 2.832196e-06
# Podle CLV platí toto jen v limitě:
a = qnorm(0.025,mean=p_chlapec,sd=sqrt(vyberovy_rozptyl/n))
b = qnorm(0.975,mean=p_chlapec,sd=sqrt(vyberovy_rozptyl/n))

cat("Výběrový průměr", vyberovy_prumer, "nespadá do intervalu od", a, "do", b, "proto nulovou hypotézu zamítneme.\n")
## Výběrový průměr 0.5139382 nespadá do intervalu od 0.5118064 do 0.5125838 proto nulovou hypotézu zamítneme.
# Nakreslíme graf N(p_chlapec, sqrt(vyberovy_rozptyl/n)).
# Červené vertikální přímky vyznačují kritické hodnoty.
# Modrý bod ukazuje hodnotu výběrového rozptylu.
x = seq(0.51, 0.515, length=100)
y = dnorm(x, mean=p_chlapec, sd=sqrt(vyberovy_rozptyl/n))
plot(x, y, type='l')
abline(v=a, col="red")
abline(v=b, col="red")
points(vyberovy_prumer, 0, col = "blue")

Zamítneme nulovou hypotézu a příjmeme alternativní. Odvodíme, že pro Českou republiku udávaný poměr narození chlapců a dívek by potřeboval upřesnit.

Pro zajímavost, ve kterém roce byla největší pravděpodobnost pro narození chlapce a ve kterém nejmenší?

rozdily = n*0 # vektor počtu narozených (chlapců-dívek) pro každý rok
divky = n*0 # vektor počtu narozených dívek pro každý rok

for (i in 1:n)
{
  rozdily[i] = data[i*2-1] - data[i*2]
  divky[i] = data[i*2]
}

max_index = which.max(rozdily)
min_index = which.min(rozdily)

max_chlapci = chlapci[max_index]
min_chlapci = chlapci[min_index]

max_divky = divky[max_index]
min_divky = divky[min_index]

cat("Největší pravděpodobnost pro narození chlapce: ", max_chlapci/(max_chlapci + max_divky), " v roce: ", 1946+max_index, "\n")
## Největší pravděpodobnost pro narození chlapce:  0.5177242  v roce:  1948
cat("Nejmenší pravděpodobnost pro narození chlapce: ", min_chlapci/(min_chlapci + min_divky), " v roce: ", 1946+min_index, "\n")
## Nejmenší pravděpodobnost pro narození chlapce:  0.5122218  v roce:  1999
narozeni_chlapce = chlapci / (chlapci + divky)

hist(narozeni_chlapce)

Příklad 2

Daly by se počty narozených dětí v desetitisících dobře modelovat Poissonovským rozdělením?

# Provedeme test dobré shody

soucty_des = round(soucty / 10000) # získáme data v desetitisících

cat("Součty v desetitisících:", soucty_des, "\n\n")
## Součty v desetitisících: 21 20 19 19 19 18 17 17 17 16 16 14 13 13 13 13 15 15 15 14 14 14 14 15 15 16 18 19 19 19 18 18 17 15 14 14 14 14 14 13 13 13 13 13 13 12 12 11 10 9 9 9 9 9 9 9 9 10 10 11 11 12 12 12 11 11 11 11 11 11 11 11
hist(soucty_des)

prihradky = 15 # hodnoty >= 15 uložíme na 15. index
lambda = mean(soucty_des) # ekvivalentně lambda = sum(soucty_des)/n
cat("lambda:", lambda, "\n")
## lambda: 13.69444
var = var(soucty_des)
cat("var:", var, "\n")
## var: 10.18701
body = 0:prihradky
p = dpois(body,lambda)
p[prihradky+1] = 1-ppois(prihradky-1,lambda)

cat("Kontrolní součet pravděpodobností:", sum(p), "\n\n")
## Kontrolní součet pravděpodobností: 1
freq = body*0

for (i in body)
{ 
  freq[i+1] = sum(soucty_des==i) # započítáme frekvence narození
}

freq[prihradky+1] = sum(soucty_des>=prihradky) # na poslední index započítáme všechna data >= 15

cat("Rozdělení do přihrádek:", freq, "\n")
## Rozdělení do přihrádek: 0 0 0 0 0 0 0 0 0 8 3 11 5 10 10 25
stopifnot(sum(freq) == n)

N = sum(freq)

T = sum((freq-p*N)^2/(p*N))
cat("Výpočet Pearsonovy statistiky podle vzorečku:", T, "\n")
## Výpočet Pearsonovy statistiky podle vzorečku: 16.56243
cat("Dosazení do distribuční funkce:", pchisq(T,prihradky), "\n\n")
## Dosazení do distribuční funkce: 0.6543122
chisq.test(x=freq,p=p)
## Warning in chisq.test(x = freq, p = p): Chi-squared approximation may be
## incorrect
## 
##  Chi-squared test for given probabilities
## 
## data:  freq
## X-squared = 16.562, df = 15, p-value = 0.3457

Tato data lze relativně dobře modelovat Poissonovským rozdělením.

Příklad 3

Proveďme lineární regresi pro počet narozených chlapců v závislosti na celkovém narození dětí a potom to samé pro dívky.

# Výpočet podle vzorečků
x = soucty
y = chlapci

xm = mean(x)
ym = mean(y)
a = sum((x-xm)*(y-ym))/sum((x-xm)^2)
b = ym - a*xm

# Výpočet podle knihovní funkce
relation <- lm(y ~ x)
summary(relation)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -695.74 -144.12    0.44  164.84  560.09 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.724e+02  1.246e+02   -2.99  0.00385 ** 
## x            5.168e-01  8.858e-04  583.39  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 235.9 on 70 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 3.403e+05 on 1 and 70 DF,  p-value: < 2.2e-16
plot(x, y)
#lines(x, a*x+b, col="blue")
abline(relation, col="red")

# To samé pro dívky

# Výpočet podle vzorečků
y = divky

ym = mean(y)
a = sum((x-xm)*(y-ym))/sum((x-xm)^2)
b = ym - a*xm

# Výpočet podle knihovní funkce
relation <- lm(y ~ x)
summary(relation)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -560.09 -164.84   -0.44  144.12  695.74 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.724e+02  1.246e+02    2.99  0.00385 ** 
## x           4.832e-01  8.858e-04  545.51  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 235.9 on 70 degrees of freedom
## Multiple R-squared:  0.9998, Adjusted R-squared:  0.9998 
## F-statistic: 2.976e+05 on 1 and 70 DF,  p-value: < 2.2e-16
plot(x, y)
#lines(x, a*x+b, col="blue")
abline(relation, col="red")