---
title: "cvič 12"
output:
  pdf_document: default
  html_notebook: default
---

# 1. příklad -- test. hypotézy, norm. veličina

```{r}
qnorm(0.95)
qnorm(0.975)
5-qnorm(0.975)
5+qnorm(0.975)
5-qnorm(0.95)
5+qnorm(0.95)
```

```{r}
u = 5 + 1.96/sqrt(10)
l = 5 - 1.96/sqrt(10)
l; u
pnorm(u,mean=4,sd=1/sqrt(10)) - pnorm(l,4,1/sqrt(10))
pnorm((u-4)*sqrt(10)) - pnorm((l-4)*sqrt(10))
```
```{r}
l2 = -10000
u2 = 5 + 1.64/sqrt(10)
l2; u2
pnorm(u2,mean=4,sd=1/sqrt(10)) - pnorm(l2,4,1/sqrt(10))
pnorm((u2-4)*sqrt(10)) - pnorm((l2-4)*sqrt(10))
```

```{r}
u3 = 10000
l3 = 5 - 1.64/sqrt(10)
l3; u3
pnorm(u3,mean=4,sd=1/sqrt(10)) - pnorm(l3,4,1/sqrt(10))
pnorm((u3-4)*sqrt(10)) - pnorm((l3-4)*sqrt(10))
```


```{r}
u = 5
pnorm(u,mean=4,sd=sqrt(10))
pnorm((u-4)/sqrt(10))
```

# 2. příklad -- chyba stroje

```{r}
qbinom(0.95,600,0.03)
pbinom(24:25,600,0.03)
```
```{r}
p=0.03
600*p + sqrt(600*p*(1-p))*qnorm(0.95)
```

# 3. příklad -- kostka


Vyřešíme dvěma způsoby (přímé dosazení do vzorce nebo použití funkce chisq.test). 
V obou případech samozřejmě vyjde totéž. 

```{r}
N=50
kostka = sample(6,N, replace=T)
#kostka = c(3,5,1,5,...)
kostka

obs = rep(0,6); obs
for(i in 1:6){
  obs[i] = sum(kostka==i)
}
obs
sum(obs)
stopifnot(sum(obs)==N)

T = sum((obs-N/6)^2/(N/6)); T
1-pchisq(T,5)    # P(hodnota T nebo vyšší)
qchisq(0.95,5)

chisq.test(obs)
```
Totéž stručněji. 

```{r}
kostka = sample(6, 100, replace=T)
obs=table(kostka); obs
chisq.test(obs)
```

Pro ilustraci: graf hustoty chi-square rozdělení s pěti stupni volnosti. Zelená plocha odpovídá 5 %, neboli 
souřadnice jejího začátku je kvantilová funkce toho rozdělení v 0.95, což je cca 11.

```{r}
x = seq(0,50,.1)
y = dchisq(x,5)
plot(x,y, type='l')
h=qchisq(0.95,5)
polygon(c(x[x>=h], max(x), h), c(y[x>=h], 0, 0), col="light green")
```
Pro kontrolu změříme pravděpodobnost chyby prvního druhu: 

```{r}
h = qchisq(0.95,5); h
fail = 0
for (i in 1:10^4){ 
  if (chisq.test(table(sample(6, 100, replace=T)))$statistic > h) { 
    fail = fail+1
  }
}
fail/10^4
```

# 4. příklad -- emaily


```{r}
moje_emaily11 = c(0,6,14,8,8,9,3,3,12,12,15,7,15,2,5,13,5,17,15,11,9,2,16,8,9,11,6,2,2,9)
#moje_emaily12 = c(13,14,3,8,5,4,12,22,8,4,5,3)
mean(moje_emaily11)
var(moje_emaily11)
hist(moje_emaily11)
day = 1:30
data = moje_emaily11[day %% 7!=1 & day %%7 != 0 & day!=17]
lambda = mean(data); lambda
var(data)
hist(data)

#data = moje_emaily12[day %% 7!=5 & day %%7 != 6 & day <= length(moje_emaily12)]
#data
#mean(data)
#var(data)
```


```{r}
n = 15
bins = 0:n
day = 1:30
data = moje_emaily11[day %% 7!=1 & day %%7 != 0 & day!=17]
lambda = mean(data); lambda
var(data)

p = dpois(bins,lambda)
p[n+1] = 1-ppois(n-1,lambda)
p
sum(p)

freq = bins*0
freq = rep(0,n+1)
for(i in bins){ freq[i+1] = sum(data==i) }
freq[n+1] = sum(data>=n)
freq

N = sum(freq); N
length(data)
stopifnot(length(data)==sum(freq))

T = sum((freq-p*N)^2/(p*N))
T
options(digits = 8)
1-pchisq(T,n)
qchisq(0.95,n)

chisq.test(x=freq,p=p)
```
Ilustrace toho, jak data sdružovat do menšího počtu přihrádek. 

```{r}
bin_ends = c(-Inf,8,10,12,Inf)  # hraniční body intervalů

day = 1:30
data = moje_emaily11[day %% 7!=1 & day %%7 != 0 & day!=17]
lambda = mean(data); lambda
var(data)
hist(data)

p = diff(ppois(bin_ends,lambda))
sum(p)

freq = bin_ends*0
for(i in 1:length(bin_ends)){
  freq[i] = sum(data<=bin_ends[i])
}
freq = diff(freq)
freq
N = sum(freq); N
length(data)

T = sum((freq-p*N)^2/(p*N))
T
1-pchisq(T,length(freq)-1)
qchisq(.95,length(freq)-1)

chisq.test(x=freq,p=p)

```

# 5. příklad -- regrese

V zadání byla řečena data pro x a y. Zde vidíme (v zakomentované části) i jak byla data vyrobena: k ideálnímu vzorci přičteme náhodný šum. Můžeme pak dobře sledovat, jak se spočtené řešení bude lišit od "ideálu". 

```{r}
x = c(1,5,9,10,13,16,20,30)
#y = c(6.1982, 12.9892, 23.8005, 23.8891, 30.0391, 35.7535, 49.0685, 63.1825) 
y = 2*x+4 + rnorm(length(x),0,3)
#y = 2*x+4 + rnorm(1,0,3)
plot(x,y)
```
```{r}
xm = mean(x)
ym = mean(y)
a = sum((x-xm)*(y-ym))/sum((x-xm)^2); a
cov(x,y)/var(x)
b = ym - a*xm; b
```
Červená je ta původní přímka, před přičtením šumu. 

```{r}
plot(x,y)
lines(x,a*x+b, col="blue")
lines(x,2*x+4, col="red")
```
Knihovní funkce na regresi ("linear model"). 
Umí např. i závislost na více proměnných (tak lze např. hledat aproximující polynom). 

```{r}
relation <- lm(y~x)

summary(relation)
relation$coefficients
```
```{r}
res = y-(a*x+b); res
```