---
title: "EDA"
author: "Andrzej Dąbrowski"
date: 
output: html_document
---

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

```{r}
library(dplyr)
library(ggplot2)
```


##Lista 1
```{r,l1z1}
load("cavendish.Rda")
cv_mean <-mean(cavendish$X1) 
cv_med <- median(cavendish$X1)
cv_mean
cv_med
g1 <- ggplot(cavendish,aes(nr,X1))
g2 <- g1 + geom_point() + 
  labs(y="pomiar [g/cm^3]")
g2
# g2 + geom_hline(yintercept = 5.515)
# g2 + geom_hline(yintercept = cv_mean)
# g2 + geom_hline(yintercept = cv_med)
g2 + geom_hline(yintercept = 5.515) + geom_hline(yintercept = cv_mean, color="red") + geom_hline(yintercept = cv_med,color="blue") + 
  annotate("segment", x = 15, xend = 18, y = 4.6, yend = 4.6,
  colour = "black") + 
  annotate("text", x = 23, y = 4.63, label = "tables                     5.515") +
annotate("segment", x = 15, xend = 18, y = 4.5, yend = 4.5,
  colour = "blue") + 
  annotate("text", x = 23, y = 4.53, label = "Cavendish median 5.460",colour = "blue") +
annotate("segment", x = 15, xend = 18, y = 4.4, yend = 4.4,
  colour = "red") + 
  annotate("text", x = 23, y = 4.43, label = "Cavendish mean    5.420",colour = "red")
```

```{r, l1z2}
wypadki <- data.frame(dzien=1:7, wypadki=c(133, 126, 159, 136, 113, 12, 20
))
wp_mean <-mean(wypadki$wypadki) 
wp_med <- median(wypadki$wypadki)
w1 <- ggplot(wypadki,aes(dzien,wypadki))
w2 <- w1 + geom_point()

w2
w2 + geom_hline(yintercept = wp_mean, color="red") + geom_hline(yintercept = wp_med,color="blue")
```

###Zad.8
**Rozrzut - odległość euklidesowa**  
```{r}
roz_eu <- function(dane){
  m <- mean(dane)
  roz <- sqrt(sum((dane-m)^2))
}
roz_czeb <- function(dane){
  roz <- (max(dane)-min(dane))/2
}
roz_taxi <- function(dane){
  m <- median(dane)
  roz <- sum(abs(dane-m))
}
roz_eu_n <- function(dane){
  n <- length(dane)
  rn <- roz_eu(dane)/sqrt(n)
}
roz_tx_n <- function(dane){
  n <- length(dane)
  rn <- roz_taxi(dane)/n
}
```

```{r}
eu <- round(c(roz_eu(cavendish$X1),roz_eu(wypadki$wypadki)),digits=2)
cz <- round(c(roz_czeb(cavendish$X1),roz_czeb(wypadki$wypadki)),digits=2)
tx <- round(c(roz_taxi(cavendish$X1),roz_taxi(wypadki$wypadki)),digits=2)
data.frame(eu=eu,czeb=cz,taxi=tx)
```

```{r}
eu_n <- round(c(roz_eu_n(cavendish$X1),roz_eu_n(wypadki$wypadki)),digits=2)
cz_n <- round(c(roz_czeb(cavendish$X1),roz_czeb(wypadki$wypadki)),digits=2)
tx_n <- round(c(roz_tx_n(cavendish$X1),roz_tx_n(wypadki$wypadki)),digits=2)
data.frame(eu=eu_n,czeb=cz_n,taxi=tx_n)
```
**Po usunięciu outliera 4.07**
```{r}
cav0 <- cavendish %>% filter(X1>4.5)
round(c(mean(cav0$X1),(max(cav0$X1)+min(cav0$X1))/2,median(cav0$X1)), digits=2)
round(c(roz_eu_n(cav0$X1),roz_czeb(cav0$X1),roz_tx_n(cav0$X1)),digits=2)
```

## Lista 2 

### Zad.8
```{r,srednia ucieta}
sredniaUcieta <- function(wektor,alfa){
  n <- length(wektor)
  if (alfa>=(n-1)/(2*n)) {
    warning("Alfa za duże")
  }
  if (alfa<(n-1)/(2*n)){
    r <- floor(alfa*n)
    i_max <- floor((n+1)/2)
    u <- alfa*n-r
    m <- n*(1-2*alfa)
    wektor <- sort(wektor)
    wagi <- rep(0,n)
    wagi[r+1] <- wagi[n-r] <- 1-u
    for (i in (r+2):i_max)  {wagi[i] <- wagi[n+1-i] <- 1}
    wyn <-(wektor %*%wagi)/(n*(1-2*alfa))
    unlist(wyn)
  }
}
sredniaUcietaV <- Vectorize(sredniaUcieta,vectorize.args="alfa")
```

###Dane
```{r}
newcomb <- c(-44,-2,16,19,21,22,23,24,24,25,26,27,28,29,29,30,31,33,34,40)
new_su <- sredniaUcietaV(newcomb,(0:9)*0.05)
new_su
require(ggplot2)
qplot((0:9)*0.05,new_su,geom =c("point", "path"))
```

```{r}
sredniaUcieta(newcomb,0.1)
```

```{r}
new_cav <- sredniaUcietaV(cavendish$X1,(0:9)*0.05)
new_cav
require(ggplot2)
qplot((0:9)*0.05,new_cav,geom =c("point", "path"))
sredniaUcieta(cavendish$X1,0.15)
```


##Lista 3

```{r, dane}
Mazowsze <- c(24,28,32,34,36,36,38,39,40,43)
Slask <- c(14,17,25,30,31,35,38,42)
```


```{r, kwantyle}
QuantileF <- function(wektor,alpha){
  wektor <- sort(wektor)
  n <- length(wektor)
  zakres <- ((alpha>=1/(2*n)) & (alpha<=(1-1/(2*n))))
  if (!zakres) stop(cat("Alpha = ",alpha, " poza zakresem"))
  if (zakres) {
  t <- alpha*n+1/2
  r <- floor(t);u <- t-r
  QF <- (1-u)*wektor[r]+u*wektor[r+1]
  }
  QF
}
```



```{r}
QuantileF(c(1,2,4,10),0.8)
QuantileF(Mazowsze,0.2187)
```





```{r}
MQu <- function(wektor){(QuantileF(wektor,0.75)+QuantileF(wektor,0.25))/2}
QQu <- function(wektor){(QuantileF(wektor,0.75)-QuantileF(wektor,0.25))/2}

MQu(Mazowsze)
MQu(Slask)
QQu(Mazowsze)
QQu(Slask)
```


```{r,rzut kwantylowy}
RzutKwantylowy <- function(wektor,k) {
  # --------------------------------------------------
  # Oblicza *k* -wymiarowy rzut kwantylowy wektora *wektor*
  # .................................................
  # Argumenty:
  # wektor: wektor danych
  # k     : rozmiar rzutu (k<=length(wektor))
  # .................................................
  # Wyniki:
  # wektor.wyn: rzut kwantylowe wektora
  #----------------------------------------------------
  
  n <- length(wektor)
  wektor <- sort(wektor)
  if (k>n) {
    stop(cat("Za duży rozmiar rzutu k = ",k))
  } 
  if (k == n) {wektor.wyn <- wektor}
  if (k<n){
    wektor.wyn <- rep(Inf,k)
    for (i in 1:k) {
      wektor.wyn[i] <- QuantileF(wektor,(2*i-1)/(2*k))
    }
  }
  wektor.wyn
}
```


```{r,przykładyrzut kwantylowy}
RzutKwantylowy(c(1,5,5,8),3)
RzutKwantylowy(c(8,5,1,5),3)
RzutKwantylowy(Mazowsze,8)
```


```{r,qqplot}
QqPlot <- function(x,y,wykres=TRUE){
  require(ggplot2)
  n <- length(x);r <- length(y)
  k <- min(n,r)
  xx <- RzutKwantylowy(x,k)
  yy <- RzutKwantylowy(y,k)
  if (wykres) {
    qq.plot <- qplot(xx,yy,xlab=deparse(substitute(x)), ylab = deparse(substitute(y)), geom =c("point", "path"))
    print(qq.plot)
    print(qq.plot+geom_smooth())
    print(qq.plot+stat_smooth(method = "lm", col = "red"))
  }
  list(x.qq=xx,y.qq=yy)
}

QqLineApx <- function(x,y) {
  qqp <- QqPlot(x,y,wykres=FALSE)
  a <- QQu(qqp$y.qq)/QQu(qqp$x.qq)
  b <- MQu(qqp$y.qq)-a*MQu(qqp$x.qq)
  qq.data <- data.frame(xx=qqp$x.qq,yy=qqp$y.qq)
  eq  <-  paste0("u = ", round(a,2), "*t + ", round(b,2))
# Plot
  sp <- ggplot(data=qq.data, aes(x=xx, y=yy)) + geom_point() + xlab(deparse(substitute(x))) + ylab(deparse(substitute(y)))
print(sp + geom_abline(intercept = b, slope =a, color="green")+
  ggtitle(eq))

  list(a=a,b=b)
}
```



```{r, przyklady qqplot}
QqPlot(Slask,Mazowsze)
sm1 <- QqPlot(Slask[-1],Mazowsze[-1])
summary(lm(sm1$y.qq~sm1$x.qq))

QqLineApx(Slask,Mazowsze)
QqLineApx(Slask[-1],Mazowsze[-1])
```

##Lista 3a

```{r dane1, echo=FALSE}
d1 <- data.frame(wzrost=c(160,180,170),masa=c(80,60,76))
d2 <- data.frame(dochod=c(7.5,13.8,15.1,19.8,20.9),inflacja=c(9.2,2.6,1.4,2.2,0.6),wyzsze=c(7,12,17,20,23))
rownames(d2) <- c("Polska","Portugalia","Hiszpania","Irlandia","Niemcy")
d1
d2

```
**Zad.1**  
```{r zad 1_3}
m1 <- apply(d1,2,mean)
m2 <- apply(d2,2,mean)
m1
m2
```
**Zad.2**  
```{r centrowanie}
centr <- function(dane, src){
  t(t(dane)-src)
} 
centr(d1,m1)
centr(d2,m2)
```
```{r macierz wariancji}
variance <- function(dane,src){
  n <- nrow(dane)
  x <- centr(dane,src)
  v <- t(x)%*%x
  v <- 1/n*v
  v
}
VX <- variance(d1,m1)
VX
VY <- variance(d2,m2)
VY
```
**Zad.3**
```{r standaryzacja}
ZX <- scale(d1)
ZX
ZY <- scale(d2)
ZY

```
**Zad.4**
```{r macierze korelacji}
CX <- cor(d1)
CX
CY <- cor(d2)
CY
```
**Zad.5**
```{r bezwladnosc}
trace <- function(data)sum(diag(data))
JVX <- trace(VX)
JVX
JVY <- trace(VY)
JVY

JCX <- trace(CX)
JCX
JCY <- trace(CY)
JCY
```
**Zad.10**

```{r rzut d1}

rzut1 <- function(x,u,a){
 z <- as.vector(unlist(x-a))
 lambda <- as.numeric(crossprod(z,u)/crossprod(u))
 y <- lambda*u+a
 y
}

```

```{r}
d2 <- function(x,y){
  z <- (x-y)^2
  d <- sum(z)/nrow(z)
  d
}
```

```{r rzuty d1 wynik}
ph <- rbind(rzut1(d1[1,],c(1,-1),c(170,70)),rzut1(d1[2,],c(1,-1),c(170,70)),rzut1(d1[3,],c(1,-1),c(170,70)))
ph
pk <- rbind(rzut1(d1[1,],c(1,-1),m1),rzut1(d1[2,],c(1,-1),m1),rzut1(d1[3,],c(1,-1),m1))
pk
```
```{r odleglosci}
d2(ph,d1)
d2(pk,d1)
```
**Zad.11**

```{r rzut d2}

rzut2  <- function(x,u1,u2){
 z <- as.vector(unlist(x))
 u11 <- as.numeric(crossprod(u1))
 u22 <- as.numeric(crossprod(u2))
 u12 <- as.numeric(crossprod(u1,u2))
 delta <- u11*u22-u12^2
 z1 <- as.numeric(crossprod(u1,z))
 z2 <- as.numeric(crossprod(u2,z))
 alfa <- (z1*u22-z2*u12)/delta
 beta <- (z2*u11-z1*u12)/delta
 y <- alfa*u1+beta*u2
 return(list(rzut=y,alfa=alfa,beta=beta))
}

```

```{r zad 11}
u1 <- as.vector(unlist(ZY[,1]))
u2 <- as.vector(unlist(ZY[,3]))
cbind(u1,u2)
x <- as.vector(unlist(ZY[,2]))
x
rzut2(x,u1,u2)
```


##Lista 4
```{r,wykres symetrii}
dane <- c(1,3,8,5,5)
wykresSymetrii <- function(dane){
  dane <- sort(dane)
  n <- length(dane)
  k <- floor((n+1)/2)
  punkty <- array(0,dim=c(2,k))
  for (i in 1:k) {
    punkty[1,i] <- (dane[n+1-i]-dane[i])/2
    punkty[2,i] <- (dane[n+1-i]+dane[i])/2
  }
  punkty <- 
    data.frame(v1=punkty[1,],v2=punkty[2,])
  gr <- ggplot(data=punkty, 
         aes(x = punkty[,1], y = punkty[,2])) 
  gr + 
  geom_point(size=2) +
  xlab("V1") + ylab("V2") +
  geom_line()
}
wykresSymetrii(dane)
```
```{r, zad 5}
dane <- 1:7
dane1 <- dane^2
dane2 <- 1-1/dane
wykresSymetrii(dane)
wykresSymetrii(dane1)
wykresSymetrii(dane2)
```


####Zad.6
```{r}
load("E:/MOJE/ANDRZEJ/dydaktyka/Uniwersytet/EAD/2015/planety.Rda")
wykresSymetrii(planety$d)
```


```{r}
load("E:/MOJE/ANDRZEJ/dydaktyka/Uniwersytet/EAD/2015/cavendish.Rda")
wykresSymetrii(cavendish$X1)
qplot(1:29,sort(cavendish$X1),geom =c("point", "path"))
qplot(cavendish$X1)
wykresSymetrii(cavendish$X1[-(1:2)])
```

```{r}
load("E:/MOJE/ANDRZEJ/dydaktyka/Uniwersytet/EAD/2015/CO2.Rda")
wykresSymetrii(CO2$co2)
qplot(CO2$co2)
```

####Współczynnik symetrii
```{r}
wspSymetrii <- function(alfa, dane){
  wsp <- (QuantileF(dane,1-alfa) + QuantileF(dane,alfa) -2*QuantileF(dane,0.5))/(QuantileF(dane,1-alfa) - QuantileF(dane,alfa))
return(wsp)
  }
```



```{r}

wspSymetrii(0.25,CO2$co2)
wspSymetrii(0.49,CO2$co2)
wspSymetrii(0.1,CO2$co2)
```

####Plotting function

<http://t-redactyl.io/blog/2016/03/creating-plots-in-r-using-ggplot2-part-9-function-plots.html>

```{r}
coefSymPlot <- function(dane, zakres=c(0.1,0.49)) {
  wsp <- function(x) {wspSymetrii(x,dane)}
  
  wyk <- ggplot(data.frame(x = zakres), aes(x = x)) +
          stat_function(fun = wsp) +
          geom_hline(yintercept = 0) +
          geom_hline(yintercept = 0.333, colour = "red") +
          geom_hline(yintercept = -0.333, colour = "blue") +
          geom_vline(xintercept = 0.25, colour = "darkgreen")
  print(wyk)
}
```

```{r,warning=FALSE}
coefSymPlot(planety$d)
coefSymPlot(cavendish$X1)
coefSymPlot(CO2$co2)
```

####Symetryzacja BoxCox
```{r}
BoxCoxExp <- function(dane){
  dane <- sort(dane)
  n <- length(dane)
  k <- floor((n+1)/2)
  md <- median(dane)
  bc <- vector(mode="logical",k)
  for (i in 1:k) {
    if (dane[i]==md) (bc[i] <- Inf)
    if (dane[n+1-i]==md) (bc[i] <- Inf)
    if (!((dane[i]==md) | (dane[n+1-i]==md)))  bc[i] <- md/(dane[n+1-i]-md) - dane[i]/(md-dane[i])
  }
  bc <- setdiff(bc,{Inf})
  return(list(bc=bc,ex=median(bc)))
}
```

####Zad. 17
```{r}
BoxCoxExp(planety$d)
BoxCoxExp(cavendish$X1)
BoxCoxExp(CO2$co2)
```

```{r,warning=FALSE}
coefSymPlot(log(planety$d))
coefSymPlot(sqrt(CO2$co2))
```

```{r,warning=FALSE}
wykresSymetrii(log(planety$d))
wykresSymetrii(sqrt(CO2$co2))
```

```{r}
wykresSymetrii(log(planety$d))+
  geom_abline(slope = 1,intercept = median(log(planety$d)),colour="red") +
  geom_abline(slope = -1,intercept = median(log(planety$d)),colour="red")

wykresSymetrii(sqrt(CO2$co2))+
  geom_abline(slope = 1,intercept = median(sqrt(CO2$co2)),colour="red") +
  geom_abline(slope = -1,intercept = median(sqrt(CO2$co2)),colour="red")
```

