#
# Zadanie 1.1
#

wektor = (1:100)^2
summary(factor(wektor %% 10))

#
# Zadanie 1.2
#

# radiany na stopnie
rad2dec <- function(x) x*pi/180
#liczymy
x = rad2dec(c(0,30,45,60,90))
data.frame(sin = sin(x), cos = cos(x), tan = tan(x), atan = 1/tan(x))

#
# Zadanie 1.3
#

paste(1:30,c("A", "B", "C"), sep=".")
paste(1:30,".",c("A", "B", "C"), sep="")

#
# Zadanie 1.4
#

daneO = read.table("http://www.biecek.pl/R/dane/daneO.csv", sep=";", h=T)
for (i in 1:ncol(daneO))
   print(paste(colnames(daneO)[i], "klasa:", class(daneO[,i]), "typ:", typeof(daneO[,i])))

#
# Zadanie 1.5
#

parzyste = 1:(nrow(daneO)/2)*2
parzyste = seq(2,nrow(daneO),2)
parzyste = (1:nrow(daneO)%%2) ==0
daneO[parzyste,]

#
# Zadanie 1.6
#

indeksy = (daneO$Wiek>50) & (daneO$Wezly.chlonne ==1)
daneO[indeksy,]

#
# Zadanie 1.7
#

(kolumny = colnames(daneO))
nchar(kolumny)

#
# Zadanie 1.8
#

skrajne <- function(x) {
 dlugosc = length(x)
 if (dlugosc<3) {
 		print("za kr�tki argument")
 		return()
 }
 sort(x)[c(1:3,(dlugosc-2):dlugosc)]
}

#
# Zadanie 1.9
#

skrajne <- function(x, ile=3) {
 dlugosc = length(x)
 if (dlugosc<ile) {
 		print("za kr�tki argument")
 		return()
 }
 sort(x)[c(1:ile,(dlugosc-ile+1):dlugosc)]
}

#
# Zadanie 1.10
#

poczatek = function(x,n) UseMethod("poczatek")
poczatek.default    = function(x,n) {
  elementow = min(n, length(x))
  x[1:elementow]
}
poczatek.matrix     = function(x,n) {
  elementow1 = min(n, dim(x)[1])
  elementow2 = min(n, dim(x)[2])
  x[1:elementow1, 1:elementow2]
}
poczatek.data.frame = poczatek.matrix

#
# Zadanie 1.11
#

gwiazdka <- function(ile = 1, max=100) {
 t = seq(0,max,0.001)
 r = 1+sin(t)
 p = t*ile
 x = r*cos(p)
 y = r*sin(p)
 plot(x,y,type="l")
}

gwiazdka()
gwiazdka(2.2)
gwiazdka(1/2)
gwiazdka(1/10)

#
# Zadanie 1.12
#

curve(2*x^3-x^2+3,0,1)
x = c(0.2, 0.8)
y = 2*x^3-x^2+3
points(x, y, pch=19, cex=1.5)
napis = paste("[", x, "," , y,"]", sep="")
text(x, y+0.05, napis)

#
# Zadanie 1.13
#

daneBT = read.table("http://www.biecek.pl/R/dane/daneBioTech.csv", sep=";", dec=",", header=T)

#
# Zadanie 1.14
#

write.table(daneBT[1:10,1:3], "c:\\maleDane.txt", sep="\t", dec=".")

#
# Zadanie 1.15
#

czyPalindrom <- function(x) {
 all(x==rev(x))
}

wektor = c(1:20, numeric(10), seq(2,40,2), seq(40,2,-2), numeric(10), 20:1)
czyPalindrom(wektor)

#
# Zadanie 1.16
#

localMin <- function(x) {
 indMin = NULL
 for (i in 1:ncol(x)) {
 	 if (class(x[,i])=="numeric") 
 	   indMin = c(indMin, which(x[,i]==min(x[,i])))
 }
 x[sort(unique(indMin)),]
}

localMin(daneBT)


#
# Zadanie 2.1
#

wiekPoziomy = cut(daneBT$Wiek, c(0,45.5, 55.5, 100), c("wiek <45", "45<= wiek <=55", "wiek >55"))
table(wiekPoziomy)
table(daneBT[,2], wiekPoziomy)
addmargins(table(daneBT[,2], wiekPoziomy))
ftable(daneBT[,2], daneBT[,3], wiekPoziomy)

#
# Zadanie 2.2
#

linie = unlist(read.table("http://www.biecek.pl/R/dane/imieniny.txt", sep="\n", stringsAsFactors=F))
length(linie)

#
# Zadanie 2.3
#

imiona = sapply(strsplit(linie, " "), function(x) x[[1]])
dlugosciImion = nchar(imiona)
table(dlugosciImion)
imiona[dlugosciImion == min(dlugosciImion)]
imiona[dlugosciImion == max(dlugosciImion)]

zenskieImiona = imiona[grep("a$",imiona)]
length(zenskieImiona)

pierwszaLitera = substr(imiona, 0,1)
summary(factor(pierwszaLitera))

imiona[grep("anna$",imiona)]

imionaBezPolskich = chartr("����󟿥ʌƣ�ӏ�","aesclnozzAESCLNOZZ", imiona)
sum(imionaBezPolskich != imiona)

#
# Zadanie 2.4
#

urodzeni30X = imiona[c(grep(" 30.X$", linie), grep(" 30.X ", linie))]
sort(urodzeni30X)
liczbyImienin = sapply(strsplit(linie, " "), function(x) length(x)-1)
linie[liczbyImienin==max(liczbyImienin)]

datyImienin = sapply(strsplit(linie, " "), function(x) x[-1])
liczebnosci = table(unlist(datyImienin))
names(liczebnosci[liczebnosci == max(liczebnosci)])

dniImienin      = unlist(sapply(strsplit(linie, "[ .]"), function(x) x[seq(2,length(x),2)]))
miesiaceImienin = unlist(sapply(strsplit(linie, "[ .]"), function(x) x[seq(3,length(x),2)]))
table(dniImienin)
table(miesiaceImienin)
plot(table(dniImienin))
plot(table(miesiaceImienin))

plot(table(tolower(unlist(sapply(strsplit(linie, "[ ]"), function(x) strsplit(x[1],""))))))

#
# Zadanie 2.5
#

with(daneBT,by(Wiek,WIT,summary))

#
# Zadanie 2.6
#

mat = matrix(c(1, 5, 3, 2, 0, 5, 1, 2, 1),3,3)
det(mat)
eigen(mat)

#
# Zadanie 2.7
#

outer(1:5,1:5,function(x,y) paste(x,"*",y,"=",x*y,sep=""))

#
# Zadanie 2.8
#

pdf("c:\\hist.pdf", 5, 5)
hist(daneBT)
dev.off()

#
# Zadanie 2.9
#

polaczenie = url("http://www.gpw.pl", open="r")
zawartoscUrla = readLines(polaczenie)
close(polaczenie)

nazwy = zawartoscUrla[grep("<td class=\"indeksy_nazwa\">", zawartoscUrla)]
wartosci = zawartoscUrla[grep("<td width=\"80\" class=\"indeksy_wart2\" align=\"right\">", zawartoscUrla)]
nazwy2 = sapply(strsplit(nazwy, "[<>]"), function(x) x[3])
wartosci2 = sapply(strsplit(wartosci, "[<>]"), function(x) x[3])

data.frame(nazwy = nazwy2[2:5],wartosci = wartosci2)

#
# Zadanie 2.10
#

# przy�pieszyli�my ponad 200 razy
Rprof("profiler.out", interval = 0.01, memory.profiling=TRUE)
N    = 10000
ind1 = 1:20
ind2 = 21:40
wyniki  = numeric(N)
macierz = .Internal(matrix(rnorm(N*40,N,40),N,40, T))
for (i in 1:N) {
	  wyniki[i] = .Internal(cor(macierz[i, ind1], macierz[i, ind2], 1, F))
}
Rprof(NULL)
summaryRprof("profiler.out",memory="both")

#
# Zadanie 2.11
#

test <- function(x) {
 wynik = try(x, silent = T)
 if (class(wynik)=="try-error") 
        cat("KLAPA\n") else cat("SUKCES\n")
}

test(cor(1))

#
# Zadanie 2.12
#

library(polynom)
wiel1 = polynomial(c(3, -1, -3, 1))
wiel2 = polynomial(c(12, -1, -12, 1))
solve(wiel1)
GCD(wiel1, wiel2)

#
# Zadanie 2.13
#

by(daneBT, daneBT[,2], head, 3)

#
# Zadanie 2.14
#

diag2 <- function(mat, n=0) {
  wymiar = min(dim(mat))-abs(n)
  if (wymiar<0) return(NULL)
  if (n==0) 
     return(diag(mat))
  if (n<0) 
     return(diag(mat[-1:n,,drop=F]))
  diag(mat[,-1:-n,drop=F])
}

diag2(matrix(1:12,3,4),1)

#
# Zadanie 2.15
#

top10 <- function() {
 wartosci = sapply(ls(.GlobalEnv), object.size)
 kolejnosc = order(as.numeric(wartosci), decreasing =T)
 wartosci[kolejnosc[1:10]]
}

top10()

#
# Zadanie 2.16
#

library("RODBC")
plikXLS <- odbcConnectExcel("c://daneBioTech.xls")
sqlTables(plikXLS)
ramka1 = sqlQuery(plikXLS, "select * from \"daneBioTech0$\"")
ramka2 = sqlQuery(plikXLS, "select * from \"daneBioTech1$\"")

#
# Zadanie 3.1
#

momenty <- function(x) {
  require(e1071)
  c(srednia = mean(x), wariancja = var(x), skosnosc = skewness(x), kurtoza = kurtosis(x))
}

momenty(rnorm(20))

#
# Zadanie 3.2
#

hist(daneO$VEGF,30)
daneO$noweVEGF = log(daneO$VEGF)
hist(daneO$noweVEGF,30, probability =T)
gest = density(daneO$noweVEGF)
lines(gest$x, gest$y, col="red",lwd=2)

#
# Zadanie 3.3
#

boxplot(daneO$Wiek~daneO$Receptory.estrogenowe)
library(vioplot)
with(daneO, vioplot(Wiek[Niepowodzenia=="wznowa"], Wiek[Niepowodzenia!="wznowa"]))

#
# Zadanie 3.4
#

attach(daneO)
library(car)
sp(Wiek, VEGF, smooth=F)
sp(VEGF, Wiek, smooth=F)
sp(Wiek, VEGF, smooth=F, log="y")
sp(Wiek,VEGF,groups=Niepowodzenia, smooth=F, log="y")

#
# Zadanie 3.5
#

x = rexp(20)
plot(ecdf(x))
fitdistr(x, "exponential")

x = rexp(200)
plot(ecdf(x))

x = rexp(2000)
plot(ecdf(x))

#
# Zadanie 3.6
#

(model1 = summary(lm(Wiek~Niepowodzenia, data = daneO)))
(model2 = summary(lm(Okres.bez.wznowy~Niepowodzenia, data = daneO)))
(model3 = summary(lm(noweVEGF~Niepowodzenia, data = daneO)))
c(Wiek = model1$coef[2,4], Okres.bez.wznowy = model2$coef[2,4], noweVEGF = model3$coef[2,4])

model <- aov(noweVEGF~Receptory.progesteronowe, data = daneO)
TukeyHSD(model)
plot(TukeyHSD(model))

#
# Zadanie 3.7
#

summary(lm(cisnienie.skurczowe~wiek, data = daneSoc))

#
# Zadanie 3.8
#

summary(lm(Kreatynina.7~Kreatynina.1+Wiek, data=daneBT))
sp(Kreatynina.7~Kreatynina.1, data=daneBT)
sp(Kreatynina.7~Wiek, data=daneBT)

#
# Zadanie 3.9
#

library(nortest)
sf.test(daneO$Wiek)
sf.test(daneO$Okres.bez.wznowy)
sf.test(daneO$VEGF)

#
# Zadanie 3.10
#

N = 1000
n = 50
pwart = numeric(N)
for (i in 1:N) {
	pwart[i] = sf.test(runif(n))$p.value
}   

alfy = seq(0,1,0.01)
moce = sapply(alfy, function(x) mean(pwart<x))
plot(alfy, moce, type="l")

#
# Zadanie 3.11
#

cor.test(daneO$Wiek,daneO$Okres.bez.wznowy)
cor.test(daneO$Wiek,daneO$VEGF)
cor.test(daneO$Okres.bez.wznowy,daneO$VEGF)

#
# Zadanie 3.12
#

nowaCena = cut(mieszkania$cena, c(0,175000,Inf), c("niska","wysoka"))
nowaPowierzchnia = cut(mieszkania$powierzchnia, c(0,45,Inf), c("mala","duza"))

library(irr)
kappa2(cbind(nowaCena, nowaPowierzchnia))

library(CoCo)
CoCoObject <- makeCoCo() 
enterTable(table(nowaCena, nowaPowierzchnia), object=CoCoObject) 
partialAssociations (": nowaCena : nowaPowierzchnia", object=CoCoObject) 

#
# Zadanie 3.13
#

inteligentnaTransformacja <- function(x) {
  require(nortest)
  wynik = list(D = lillie.test(x)$statistic, trans="bez zmian", x = x)
  nx = log(x)
  if (lillie.test(nx)$statistic < wynik$D) {
     wynik = list(D = lillie.test(nx)$statistic, trans="log", x = nx)
  }
  nx = sqrt(x)
  if (lillie.test(nx)$statistic < wynik$D) {
     wynik = list(D = lillie.test(nx)$statistic, trans="sqrt", x = nx)
  }
  nx = 1/x
  if (lillie.test(nx)$statistic < wynik$D) {
     wynik = list(D = lillie.test(nx)$statistic, trans="ivers", x = nx)
  }
  wynik
}

inteligentnaTransformacja(rnorm(100,5))
inteligentnaTransformacja(1/rnorm(100,5))
inteligentnaTransformacja(rnorm(100,5)^2)
inteligentnaTransformacja(exp(rnorm(100,5)))

#
# Zadanie 3.14
#

(kontr = cbind(c(3, -1, -1, -1), c(1, 1, -1, -1), c(1, 1, 1, -3)))
summary(lm(noweVEGF~Receptory.estrogenowe, data=daneO, contrasts= list(Receptory.estrogenowe=kontr)))

#
# Zadanie 3.15
#

summary(lm(daneBT[,6]~daneBT[,2]+daneBT[,3]))
summary(lm(daneBT[,6]~daneBT[,2]*daneBT[,3]))
summary(lm(daneBT[,6]~daneBT[,2]:daneBT[,3]))

interaction.plot(daneBT[,3], daneBT[,2], daneBT[,6])

#
# Zadanie 3.16
#

bezpiecznyTest <- function(x1, x2, alpha) {
  require(nortest)
  pwart = sf.test(c(scale(x1),scale(x2)))$p.value
  if (pwart<alpha) {
   return(wilcox.test(x1,x2))
  } 
  t.test(x1,x2)
}

# blad pierwszego rodzaju

N = 100000
n = 10
alfa = 0.1
pwart = numeric(N)
for (i in 1:N) {
	pwart[i] = bezpiecznyTest(rchisq(n,2), rchisq(n,2), alfa)$p.value
}
mean(pwart<alfa)

# moc brak normalnosci

N = 10000
n = 10
alfa = 0.1
pwart1 = pwart2 = pwart3 = numeric(N)
for (i in 1:N) {
	x1 = rchisq(n,2)
	x2 = rchisq(n,2)+1
	pwart1[i] = bezpiecznyTest(x1,x2, alfa)$p.value
	pwart2[i] = t.test(x1,x2)$p.value
	pwart3[i] = wilcox.test(x1,x2)$p.value
}
mean(pwart1<alfa)
mean(pwart2<alfa)
mean(pwart3<alfa)

# moc, rozklady normalne

N = 10000
n = 20
alfa = 0.1
pwart1 = pwart2 = pwart3 = numeric(N)
for (i in 1:N) {
	x1 = rnorm(n,2)
	x2 = rnorm(n,2)+0.7
	pwart1[i] = bezpiecznyTest(x1,x2, alfa)$p.value
	pwart2[i] = t.test(x1,x2)$p.value
	pwart3[i] = wilcox.test(x1,x2)$p.value
}
mean(pwart1<alfa)
mean(pwart2<alfa)
mean(pwart3<alfa)

#
# Zadanie 4.1
#

library(gplots)
attach(daneSoc)
hist2d(cisnienie.skurczowe, cisnienie.rozkurczowe, nbins=15)

#
# Zadanie 4.2
#

attach(daneO)
sunflowerplot(daneO$Wiek[!is.na(daneO$Nowotwor)], na.omit(daneO$Nowotwor))

#
# Zadanie 4.3
#

pairs(mieszkania)
library(YaleToolkit)
gpairs(mieszkania)

#
# Zadanie 4.4
#

library(ks)
danetmp = daneSoc[, c(6,7)]
mcov = Hscv(danetmp)
gest2d = kde(danetmp, H=mcov)
plot(gest2d, cont=c(5,25,50,75,95))

#
# Zadanie 4.5
#

macierz = cor(daneBT[,c(1,4:13)])
library(ellipse)
plotcorr(macierz)
heatmap(macierz, symm=T)

#
# Zadanie 4.6
#

matplot(c(1,3,7),t(daneBT[1:5,4:6]), type="l")

#
# Zadanie 4.7
#

dwaWykresy <- function(x1, x2, col1="black", col2="red" , ...) {
 x = 1:length(x1)
 plot(x1, type="l", col=col1, ...)
 szerokosc = diff(range(x2))
 par(usr=c(0,length(x2),min(x2)-0.04*szerokosc,min(x2)+1.04*szerokosc))
 lines(x2, col=col2, ...)
 axis(side=4)
}

x0 = rnorm(500)
x1 = cumsum(x0)
x2 = cumsum(rnorm(500)-x0/4)*100

dwaWykresy(x1,x2, lwd=2)


#
# Zadanie 4.8
#

par(mfrow=c(1,2), xpd=NA)
x = runif(10)
y = runif(10)

plot(x[1:5],y[1:5], pch=8,lwd=2,cex=1.5,xlab="",ylab="",xlim=c(0,1),ylim=c(0,1))
plot(x[6:10],y[6:10], pch=1,lwd=2,cex=1.5,xlab="",ylab="",xlim=c(0,1),ylim=c(0,1))

par(usr=c(-0.04,1.04,-0.04,1.04))

margines = 2*par()$plt[1]/diff(par()$plt[1:2])

arrows(x[6:10],y[6:10],x[1:5]-1-margines,y[1:5], code=3, angle=15, col="red")