---
title: "Zestawy"
author: "Author name"
date: "6 czerwca 2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown (PRZKLAD)

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r aid}
summary(Aids2)
```

## R Markdown (END OF PRZKLAD)


## Zestaw 1 zadanie 12
Wykonaj histogram i wykres pudełkowy dla danych: aid (miesięczne
płatności na program federalny w USA), crime (współczynnik
przestępczości dla stanów w USA (lata 1983, 1993)) oraz south
(współczynnik zabójstw w amerykańskich miastach) z pakietu UsingR.
Który z tych zbiorów danych jest asymetryczny, a który symetryczny? W
jakim zbiorze istnieją obserwacje odstające?

## Zestaw 2 zadanie 16
Niech X_1,...,X_n - będą zmiennymi losowymi o rozkładzie jednostajnym w przedziale [0, 1].
Oszacować prawdopodobieństwo, że wartość średnia z próbki n- elementowej jest mniejsza niż mediana.
Obliczenia wykonać dla losowo wybranych k=100000 próbek dla n=30. Wynik zapisać w postaci tabeli.

# Wyznaczanie przedziaіu ufnoњci dla odchylenia standardowego

f3=function(x,alpha){
  n=length(x)
  chi1=qchisq(1-alpha/2,n-1)
  chi2=qchisq(alpha/2,n-1)
  x_l=n*var(x)/chi1
  x_p=n*var(x)/chi2
  list(lewy_koniec=x_l,prawy_koniec=x_p,lewy=sqrt(x_l),
       prawy=sqrt(x_p))
}

# Wyznaczanie przedziaіu ufnoњci dla proporcji

f4=function(n,m,alpha){
  p=m/n
  u_alpha=qnorm(1-alpha/2)
  x_l=p-u_alpha*sqrt(p*(1-p)/n)
  x_p=p+u_alpha*sqrt(p*(1-p)/n)
  list(lewy=x_l,prawy=x_p)
}
f4(500,100,0.05)
f4(5000,1000,0.05)
f4(50000,10000,0.05)

# Wykonamy symulacjк wyjaњniaj№c№ czym jest 95% przedziaі ufnoњci dla wartoњci sredniej.
# W pierwszej kolejnoњci wygenerujemy dwie populacje licz№ce po 6000 elementуw.
# pierwsz№ x1 - z rozkіadu N(2,1) a drug№ x2 - z rozkіadu jednostajnego w przedziale
# [-1,1].
x1<-rnorm(6000,2,1)
x2<-runif(6000,-1,1)
# Nastкpnie z tych zbiorуw bкdziemmy 100 razy losowaж prуbki po 200 elementуw.
# Dla kaїdej wylosowanej 200 elementowej prуbki wyznaczymy przedziaі ufnoњci
# i ich koсce zapiszemy w wektorach l i r
# Dobr№ zasad№ jest na wstepie utworzenie wektorуw l i r 200 elementowych
# zawieraj№cych NA. Pozwala to kontrolowaж w trakcie symulacji kroki ktуre juї zostaіy
# wykonane
k=200 # liczebnoњж pojedynczej prуbki
n<-100 # liczba iteracji
l<-rep(NA,n) # generowanie wektora gdzie bкd№ przechowywane lewe koсce przedziaіow ufnoњci
r<-rep(NA,n) # generowanie wektora gdzie bкd№ przechowywane prawe koсce przedziaіow ufnoњci

# Obliczenia dla zbioru x1.

for(i in 1:100)
{
  print(i)
  a1<-sample(1:6000,k,replace=F) # generowanie k indeksуw dla i-tej prуbki
  l[i]<-t.test(x1[a1],alternative="two.sided",mu = 0,conf.level=0.95)$conf.int[1]
  r[i]<-t.test(x1[a1],alternative="two.sided",mu = 0,conf.level=0.95)$conf.int[2]
}
# Okreњlimy ktуre ze 100 przedziaіуw ufnoњci zawiera wrtoњж њredni№ populacji x1
mean((l<=mean(x1,na.rm=T))&(mean(x1,na.rm=T)<=r))
# zobrazujemy to graficznie
win.graph()
plot(c(l,r),c(1:100,1:100),type="n",xlab="Przedziaі ufnoњci",ylab="Prуbki",
main="Prуbki losowane ze zbioru o rozkіadzie N(2,1)")
s<-seq(1,100,length=100)
for (i in 1:100)
{ifelse((l[s][i]<=mean(x1,na.rm=T))&(mean(x1,na.rm=T)<=r[s][i]),b<-1,b<-2)
  lines(c(l[s][i],r[s][i]),c(s[i],s[i]),col=b)}
lines(rep(mean(x1,na.rm=T),100),1:100)

# Obliczenia dla zbioru x2.

for(i in 1:100)
{
  print(i)
  a1<-sample(1:6000,k,replace=F) # generowanie k indeksуw dla i-tej prуbki
  l[i]<-t.test(x2[a1],alternative="two.sided",mu = 0,conf.level=0.95)$conf.int[1]
  r[i]<-t.test(x2[a1],alternative="two.sided",mu = 0,conf.level=0.95)$conf.int[2]
}
#
# Okreњlimy ktуre ze 100 przedziaіуw ufnoњci zawiera wrtoњж њredni№ populacji x2
mean((l<=mean(x2,na.rm=T))&(mean(x2,na.rm=T)<=r))
# zobrazujemy to graficznie
win.graph()
plot(c(l,r),c(1:100,1:100),type="n",xlab="Przedziaі ufnoњci",ylab="Prуbki",main=
       "Prуbki losowane ze zbioru o rozkіadzie jednostajnym w [-1, 1]")
s<-seq(1,100,length=100)
for (i in 1:100)
{ifelse((l[s][i]<=mean(x2,na.rm=T))&(mean(x2,na.rm=T)<=r[s][i]),b<-1,b<-2)
  lines(c(l[s][i],r[s][i]),c(s[i],s[i]),col=b)}
lines(rep(mean(x2,na.rm=T),100),1:100)
#
# Symulacja przedziaіуw ufnoњci dla proporcji z wykorzystaniem funkcji
# matplot z pakietu graphics.
# Piкdziesi№t razy rzucamy 20 razy monet№.
library(graphics)
m = 50; n=20; p=0.5;alpha=0.1 # p- prawdopodobieсstwo wyrzucenia orіa
ls=rbinom(m,n,p) # wektor ktуry zawiera liczbк orіуw z 50 doњwiadczeс.
matplot(rbind(f4(n,ls,alpha)[[1]],f4(n,ls,alpha)[[2]]),rbind(1:m,1:m),type="l",lty=1)
abline(v=p)

## Zestaw 3 zadanie 4
Czy istnieje różnica między grupami dotyczącymi końcowego pomiaru i jaki jest rozmiar efektu?
(Wskazówka: użyj jednokierunkowej ANOVA)


Testowanie przeprowadza siк za pomoc№ funkcji **anova()** porуwnuj№c model ktуry uwzglednia wszystkie zmienne niezaleїne z modelami ktуre zawieraj№ mniej zmiennych niezaleїnych.
Dla regrsji liniowej jednej zmiennej model liniowy porуwnuje siк z modelem ktуry uwzglкdnia tylko staі№.

Okreњlmy model zaleїny tylko od staіej:

```{r}
eruptions0.lm = lm(eruptions ~ 1, data=faithful)
```

Porуwnamy za pomoc№ funkcji **anova()** modele eruptions0.lm i eruptions.lm:

```{r}
anova(eruptions0.lm,eruptions.lm)
```

Zauwaїmy, їe wartoњж statystyki **F** i **p-value** wyznaczone powyїej s№ takie same jak wyznaczone za pomoc№ funkcji **summary(eruptions)**.
Poniewaї wartoњж **p-value** jest duїo mniejsza od 0.05 wiкc model **eruptions.lm** jest lepszy niї model **eruptions0.lm**.


## Zestaw 4 zadanie 12
Obserwując liczbę awarii w sieci wodno-kanalizacyjnej w ciągu 100 dni w pewnym rejonie miasta otrzymano dane:
Dzienna liczba awarii
0
1
2
3
4
Liczba dni
13
32
27
18
10
a) Na poziomie ufności 1 - a =0,9 oszacować metodą przedziałową średnią dzienną liczbę awarii w l losowo
wybranym dniu.
b) Na poziomie ufności 1 - a =0,95 oszacować metodą przedziałową wariancję dziennej liczby awarii w sieci
wodno-kanalizacyjnej.
c) Na poziomie istotności a = 0,05 zweryfikować hipotezę, że średnia dzienna liczba awarii w sieci wodno-kanalizacyjnej
jest równa 1,5.


# Przykіad 6.
#
# Administrator sieci komputerowej zebraі dane na temat liczby awarii sieci podczas ostatnich
# 500 dni. Uzyskaі nastкpuj№ce wyniki:
#  l_awarii   l_dni
#   0          119
#   1          182
#   2          123
#   3           48
#   4           17
#   5            6
#   6            5

# Czy na poziomie istotnoњci 0.05 moїna s№dziж, їe liczba awarii ma rozkіad Poissona?
# Niech X oznacza liczbк awarii. Mamy nastкpuj№ce hipotezy:
#  H0 : X ~ pois(lambda) - zmienna losowa X ma rozkіad zgodny z rozkіadem Poissona
# H1 : ~ (X ~ pois(??) ).
# Nie znamy parametru lambda w rozkіadzie Poissona, ale wiemy, їe moїna go oszacowaж przez њredni№
# (jako wartoњж oczekiwana) :
l_awarii=0:6  # wektor okreњlaj№cy liczbк awarii
l_dni=c(119, 182, 123, 48, 17, 6, 5) # wektor okreњlaj№cy liczbк dni poszczegуlnych awarii
dane=rep(l_awarii, l_dni) # tworzymy wektor danych
lambda=mean(dane)  # wyliczamy wartoњж њredni№
prawd=c(dpois(0:5, lambda), ppois(5, lambda, lower.tail=F)) # obliczamy prawdopodobieсstwa
# dla poszczegуlnych liczby awarii
chisq.test(table(dane), p=prawd)
#

## Zestaw 5 zadanie 10



Zadanie 10. a. Wczytaj zbiór danych star z pakietu faraway i zapoznać się ze zmiennymi tego zbioru. b.
Zbadać zależność zmiennej light od zmiennej temp zarówno liniową jak i wielomianową
drugiego stopnia. c. Wykonać stosowne wykresy i analizę dopasowania.



#### Miara dopasowania

Najbardziej popularn№ miar№ dopasowania jest wspуіczynnik determinacji. Jest to liczba z przedziaіu $[0,\,1]$ i oznacza siк j№ $R^2$, wartoњж $R^2=1$  oznacza doskonaіe dopasowanie, natomiast wartoњж $R^2=0$ oznacza brak powi№zania liniowego miкdzy zmiennymi. Zanim okreslimy wielkoњж $R^2$ wprowadџmy nastкpuj№ce oznaczenia:
$$\sum_{i=1}^n(y_i-\bar{y})^2 -\quad  caіkowita\,\, suma\,\, kwadratуw\,\, (CSK),$$
$$\sum_{i=1}^n(\hat{y_i}-\bar{y})^2 -\quad wyjaњniona\,\, suma\,\, kwadratуw\,\, przez\,\,model\,\, (WSK),$$
$$\sum_{i=1}^ne_i^2 \, -\quad  resztowa\,\, suma\,\,kwadratуw\,\,(RSK).$$
Miкdzy tymi wielkoњciami zachodzi nastкpuj№ca zaleїnoњж:
$$\sum_{i=1}^n(y_i-\bar{y})^2=\sum_{i=1}^n(\hat{y_i}-\bar{y})^2 + \sum_{i=1}^ne_i^2.$$
Jako wspуіczynnik determinacji przyjmujemy stosunekk zmiennoњci wyjaњnionej do zmiennoњci caіkowitej.
Zatem
$$R^2=\dfrac{WSK}{CSK}.$$
Wspуіczynnik $R^2$ mierzy nam, jaka czкњж ogуlnej zmiennoњci zmiennej zaleїnej jest wyjaњniona przez regresjк liniow№.
Im wiкksze $R^2$, tym lepiej. Nie naleїy jednak przesadzaж. Doі№czenie bowiem nowej zmiennej do istniej№cego modelu zawsze powoduje zwiкkszenie $R^2$. Naszym celem nie jest uzyskanie jakk najwiкkszego $R^2$, a znalezienie zwi№zku miкdzy X i Y z rzetelnymi ocenami parametrуw.

#### Klasyfikacja obserwacji nietypowych

Poprawne wnioskowanie i poprawna diagnoza dopasowanego modelu jest niemoїliwa bez wykrycia i analizy wpіywu na model tzw. obserwacji nietypowych.
Obserwacje takie mog№, choж nie musz№, wywieraж istotny wpіyw na wyniki.
Wњrуd takich obserwacji bкdziemy wyrуїniaж: punkty **wysokiej dџwigni**, **punkty odstaj№ce** (oddalone), **punkty wpіywowe**.

Wielkoњж $h_i$ okreњlon№ nastкpuj№co:
$$h_i=\frac{1}{n}+\frac{x_i-\bar x}{\sum_j^n(x_j-\bar x)^2}$$
bкdziemy nazywaж **dџwigni№** $i$-tej obserwacji.

Dџwignia okreњla wariancjк i-tej reszty (wiкksza dџwignia powoduje mniejsz№ wariancjк). Wiкksz№ dџwigniк maj№ obserwacje,
ktуre znajduj№ siк na skraju zakresu zmiany zmiennej niezaleїnej, natomiast wartoњж zmiennej zaleїnej zupeіnie nie wpіywa na wielkoњж dџwigni. Najmniejsz№
dџwigniк (rуwn№ $\frac{1}{n}$) ma obserwacja, ktуrej wartoњж $x_i$ pokrywa siк z wartosci№ њredni№ $\bar x$. Kres gуrny dla wartoњci dџwigni jest rуwny jednoњci.

Zachodzi nastepuj№ca zaleїnoњж
$$Var(e_i)=(1-h_i)Var(y_i).$$
Obserwacjк o numerze $i, \,i = 1, . . . , n$, nazywamy **obserwacj№ odstaj№c№**, jeњli
$$\frac{\vert e_i\vert}{\sqrt{Var(y_i)\cdot (1-h_i)}}>2.$$
**Odlegіoњci№ Cooka** dla obserwacji o numerze $i,\, i = 1, . . . , n,$ okreњlamy nastepuj№co:
$$D_i =\frac{e_i^2h_i}{2Var(y_i)\cdot (1-h_i)^2}.$$
Obserwacjк nazywamy **wpіywow№**, jeњli dla niej $D_i > 1$.

Podczas analizy modelu liniowego interesuje nas, czy zaleїnoњж pomiкdzy zmiennymi X i Y jest istotna. Formuіuj№c to pytanie w jкzyku modeli liniowych, otrzymujemy hipotezк zerow№ dotycz№c№ wspуіczynnika $\beta_1$:
$$ H_0:\,\,\beta_1\,=\,0$$
oraz hipotezк alternatywn№
$$ H_1:\,\,\beta_1\,\neq\,0.$$
Jeїeli w wyniku testowania za pomoc№ **t-testu** odrzucimy hipotezк zerow№ to oznacza, їe pomiedzy zmiennymi X i Y jest istotna zaleznoњж.
Oprуcz testu **t** stosuje siк globalny test **F (Fiszera-Sendecora)**. Test ten weryfikuje nam rуwnoczeњnie trzy rуwnowaїne hipotezy zerowe:
$$ H_0:\,\,\beta_1\,=\,0 \,-\quad  istotnoњж\,\,wspуіczynnika\,\,kierunkowego$$
$$ H_0:\,\,R^2\,=\,0\,- \quad istotnoњж\,\,wspуіczynnika\,\,R^2,$$
$$ H_0:\,\,y=\beta_1x+\beta_0\,-\quad istotnoњж\,\,liniowego\,\,zwi№zku.$$
Przy weryfikacji **F testem** wykorzystuje siк nastepuj№c№ statystykк
$$F=\dfrac{WSK \cdot(n-2)}{RSK}.$$



Poniїsze polecenie pozwala przedstawiж dane wraz z dopasowan№ lini№ regresji:

```{r}
attach(faithful)
plot(waiting,eruptions)
abline(eruptions.lm)
```



Z poniїej przedstawionymi modelami jest coњ nie tak.

```{r}
mod2=lm(y2~x2,anscombe)
plot(y2~x2,anscombe)
abline(mod2)
summary(mod1)
print(plot(mod2, which=1:6))
```

W tym przypadku ewidentnie lepszym rozwiazaniem byіby model kwadratowy.

```{r}
mod3=lm(y3~x3,anscombe)
plot(y3~x3,anscombe)
abline(mod3)
summary(mod3)
print(plot(mod3, which=1:6))
```

W tym przypadku jedna obserwacja jest odstajaca i zaburza oceny modelu.

```{r}
mod4=lm(y4~x4,anscombe)
plot(y4~x4,anscombe)
abline(mod4)
summary(mod4)
print(plot(mod4, which=1:6))
```

W tym przypadku trend liniowy zostaі wygenerowany przez jedn№ obserwacjк i nie ma nic wspуlnego z liniow№ zaleїnoњci№ zmiennych.

#### Rysowanie krzywej regresji i przedziaіu ufnoњci dla wartoњci przewidywanych przez model

```{r}
library(faraway)
data(corrosion)
head(corrosion)
attach(corrosion)
model3=lm(loss~Fe,data=corrosion)
grid=seq(0,3,by=0.1)
p=predict.lm(model3,data.frame(Fe=grid),se=TRUE)
q=qt(0.975,11)
matplot(grid,cbind(p$fit, p$fit-q*p$se, p$fit+q*p$se), type="l",lty=c(1,2,2),
        xlab="Iron Content", ylab="loss")
rug(Fe)
points(Fe,loss)
```