###################################################################
# Tutorial de Introdução ao R
###################################################################

O programa R é gratuito e de código aberto que propicia excelente ambiente para 
análises estatísticas e com recursos gráficos de alta qualidade. Detalhes sobre 
o projeto, colaboradores, documentação e diversas outras informações podem ser 
encontradas na página oficial do projeto em:             
http://www.r-project.org. 

O programa pode ser baixado livremente pela internet. Há alguns espelhos 
(mirrors) brasileiros da área de downloads do programa chamada de CRAN 
(Compreensive R Arquive Network), que pode ser acessado em 
http://cran.br-r-project.org 

Além dos materiais disponíveis na página do programa há também Tutoriais de 
Introdução ao R disponível em:
http://leg.ufpr.br/~paulojus/embrapa/Rembrapa/Rembrapa.pdf 


O editor de texto Tinn-r também é um programa livre, e pode ser encontrado na 
página:
http://sourceforge.net/projects/tinn-r

O editor de texto RStudio também é um programa livre, e pode ser encontrado na 
página:
https://www.rstudio.com/products/rstudio/download/#download


Iniciar o R
O símbolo > indica a linha de comando

help.start()

?mean
help(mean)
args(mean)
?rnorm
args(rnorm)

1+2+3    #somando estes números, obtém-se a resposta marcada com [1]
2+3*4    #um pouquinho mais complexo, prioridade de operações (multiplicação primeiro)
3/2+1
4*3**3   #potências são indicadas por ** ou ^
4*3^3
sqrt(2)
sin(2)
cos(pi)
acos(-1)
exp(2)
log(2)
log10(2)
gamma(4)   
factorial(3)  
sqrt(sin(45*pi/180))

# Armazenando valores em objetos do R
x <- 10
X = 9
x
X
y <- 11
c(x,X,y)

x <- 1:100
x

y+x

a <- c(2,3,5,7,11)
b <- c(a,2,2)
a
b
a[5]

a[4] <- 6
a

length(a)
length(b)

is.vector(a)
is.matrix(a)
is.numeric(a)
is.character(a)

xx <- 1:10
xx
xy <- seq(1,10)
xy
xz <- seq(10,1)
xz
xw <- seq(1,10,2)
xw
xq <- seq(10,1,-2)
xq
xx; xy; xz; xw; xq;
xx[2]
xx[2:4]

xx[xx < 5]    # retorna valores numéricos
xx[xx >= 8]
xx[xx > 8] 
xx > 5        # retorna valores lógicos (V ou F)

rep(1,10)
rep(c(1,2),10)
rep(4:1,1:4)
rep(c(23, 32, 42), c(3, 1, 2))

ls()
remove(y)
ls()
remove(list=ls())
ls()            

x <- (1:100)^2
x
summary(x)

# Exercícios
#1) Defina a amostra a = (10,10,9,8,8,7,6,6,5,4,4,3).
#2) Encontre as estatísticas descritivas desta amostra.
#3) Ordene a amostra de forma crescente.
#4) Calcule  \sqrt{sen(a^2)+exp(a)}

#Resposta
a<-rep(10:3,rep(2:1,4)); a;
summary(a)
sort(a)
sqrt(sin(a^2)+exp(a))

###################################################################
# Segunda aula
###################################################################
x <- (1:100)^2
x
?plot

plot(x)
plot(x,type="l")
plot(x,type="b")
plot(x,type="o")
plot(x,type="s")
plot(x,type="c")
plot(x,type="h")

x <- 1:20
y <- x**3
plot(x,y)
points(rev(x),y,col="red")
lines(x,8000-y)

plot(x,y, xlab="Eixo X aqui",ylab="Eixo Y aqui", lwd=4, col="blue", type="l")

plot(x,y, xlab="Eixo X aqui",ylab="Eixo Y aqui", lwd=4, lty=5, col="blue", type="l")
title("Título vai aqui \n (e note a acentuação!!!)")
text(6,4000,"Texto em qualquer lugar")

par(mfrow=c(2,2))
plot(x,y)
plot(x,y)
plot(x,2*y)
plot(x,log(y))

par(mfrow=c(3,1))
x<-c(2,2,2,2,2,3,3,3,4,4,5,5)
hist(x)
barplot(table(c(2,2,2,2,2,3,3,3,4,4,5,5)))
barplot(table(c(2,2,2,2,2,3,3,3,4,4,5,5)), hor=T)

y<-summary(x)
y
names(y)

# ggplot2 exemplos
library(ggplot2) 

data(mtcars)
mtcars
names(mtcars)
summary(mtcars$gear)
mtcars$gear
length(mtcars$gear)

# criando fatores com valores "lables" 
mtcars$gear <- factor(mtcars$gear,levels=c(3,4,5),labels=c("3gears","4gears","5gears")) 
mtcars$gear
mtcars$am <- factor(mtcars$am,levels=c(0,1),labels=c("Automatic","Manual")) 
mtcars$am
mtcars$cyl <- factor(mtcars$cyl,levels=c(4,6,8),labels=c("4cyl","6cyl","8cyl")) 
mtcars$cyl

# Densidade Kernel para mpg agrupado por "number of gears" 
qplot(mpg, data=mtcars, geom="density", fill=gear, alpha=I(.5), 
      main="Distribution of Gas Milage", xlab="Miles Per Gallon", 
      ylab="Density")

# Scatterplot de mpg vs. hp para cada combinação de "gears" e "cylinders" 
qplot(hp, mpg, data=mtcars, shape=am, color=am, 
      facets=gear~cyl, size=I(3),
      xlab="Horsepower", ylab="Miles per Gallon") 

# Regressões separadas de "mpg" sobre "weight" para cada "number of cylinders"
qplot(wt, mpg, data=mtcars, geom=c("point", "smooth"), 
      method="lm", formula=y~x, color=cyl, 
      main="Regression of MPG on Weight", 
      xlab="Weight", ylab="Miles per Gallon")

# Boxplots de "mpg" por "number of gears" 
qplot(gear, mpg, data=mtcars, geom=c("boxplot", "jitter"), 
      fill=gear, main="Mileage by Gear Number",
      xlab="", ylab="Miles per Gallon")

# Customizando os gráficos
p <- qplot(hp, mpg, data=mtcars, shape=am, color=am,facets=gear~cyl, main="Scatterplots of MPG vs. Horsepower",
           xlab="Horsepower", ylab="Miles per Gallon")

# Fundo branco e linhas de gride pretas
p + theme_bw()

# Alterando novamente
p + theme(axis.title=element_text(face="bold.italic", size="12", color="brown"), legend.position="top")


#Exercício

1) Contrua o gráfico de Y=3*X+4, em que X varia de 0 a 10, com o 
formato de  linha e na cor azul.
2) Defina o nome dos eixos X, como "Eixo X" e Y, como "Eixo Y".
3) Insira um título.

#Resposta
X<- 0:10
Y<- 3*X+4
plot(X,Y, type="l", col="blue")
plot(X,Y, type="l", col="blue", xlab="Eixo X",ylab="Eixo Y")
title("Título do gráfico")

###################################################################
# Terceira aula
###################################################################
# Criando uma lista
person <- list(age=21,name='Fred',score=c(65,78,55))
person
person[[2]]
person$name
names(person)
person$score
person$score[1]
person$score[3]

#  Criando Matrizes
x <- 1:12
xmat <- matrix(x,ncol=3)    # automaticamente nrow = 4
xmat 
dim(xmat)
nrow(xmat)
ncol(xmat)
dimnames(xmat)
dimnames (xmat) <- list(c("L1","L2","L3","L4"),c("C1","C2","C3"))
dimnames (xmat)
xmat
xmat[c("L1","L3"),c("C1")]
mode(dimnames(xmat))

matrix(x,ncol=3,byrow=T)

m21 <- c(1:5);m21
m22 <- c(6:10);m22
m23 <- c(11:15);m23
plot(m22,m21)

m2 <- cbind(1:5,6:10,11:15) # funções que adicionam linhas ou colunas
m2
plot(m2[,1],m2[,2])

mm2 <- cbind(m21,m22,m23)
mm2

m3<-rbind(1:5,5:1)
m3
is.matrix(m3)

x1 <- matrix(1:6,ncol=2)
x2 <- matrix(6:1,ncol=3)
x1;x2;

x1 %*% x2 # multiplicação de matrizes

t(x1)*x2  # multiplicação de elementos de matrizes

x1[1,]
x1[,2] 

# Criando arranjos
ar1 <- array(1:24, dim = c(3, 4, 2)) 
ar1
ar1[1,1,1]
ar1[1,1,2]
ar1[,1:2,1]

sum(ar1[,1,1])
colSums(ar1[,,1])

###################################################################
# Quarta aula
###################################################################
# Análise de dados - voltando ao exemplo dos carros
?mtcars
data(mtcars)
mtcars
dim(mtcars)
mtcars[1:5,]
names(mtcars)
is.list(mtcars)
mtc <- mtcars[,c(1,2,4,6,9,10)]
mtc
mtc[1:5,]
names(mtc)

# Tabela de frequências absolutas
?table
tcyl <- table(mtcars$cyl)
tcyl
sum(tcyl)

cam <- table(mtcars$am)
cam

# Gráfico de frequências absolutas
?barplot
barplot(tcyl)

# Gráfico de tortas
?pie
pie(tcyl)

# Tabela de frequências relativas
table(mtcars$cyl)
tcyl_r <- 100* (tcyl)/length(mtcars$cyl)
tcyl_r
sum(tcyl_r)

# Gráfico de frequências relativas
barplot(tcyl_r)

# Agora dados contínuos
mpg <- mtcars$mpg
mpg
range(mpg)
length(mpg)
summary(mpg)
table(cut(mpg, br=seq(10,35,5))) # tabela de frequência absoluta com intervalos a cada 5 mpg

# Histograma
?hist
hist(mpg) # Histograma com frequencias absolutas = contagem
hist(mpg, freq=FALSE, ylim=c(0,0.1)) # Histograma com frequencias relativas/intervalos

?boxplot
boxplot(mpg, col="bisque", main = "consumo em galões por milhas")
summary(mpg)

## Descrição bivariada
?mtcars
names (mtcars)
am <- mtcars$am
am
cyl <- mtcars$cyl
cyl

table(am, cyl)
barplot(table(am, cyl), legend.text = c("automatico","manual"), main = "gráfico bidimensional Cyl vs Gear",legend = rownames(mtcars$am))

## Como obter medidas resumo do rendimento para cada tipo de câmbio
hp <- mtcars$hp;hp
tapply(mpg, am, summary)
tapply(hp, am, summary)

plot(hp, am)
plot(mpg, am)
plot(mpg, hp)

###################################################################
# Quinta aula
###################################################################
# Distribuições de probabilidade
x <- runif(100) # 100 números aleatórios

x<-rnorm(1000,2,3)
hist(x, freq=F)
curve(dnorm(x,2,3), col=2, lty=1, lwd=2, add=T)
density(Y)

# rnorm(n,mean=0,sd=1) Gera n números
# dnorm(q,mean=0,sd=1) Densidade para o q-ésimo quantil
# pnorm(q,mean=0,sd=1) Probabilidade para o q-ésimo quantil
# qnorm(p,mean=0,sd=1) Quantil para a probabilidade p

X<-rnorm(100,0,1)
summary(X)
plot(X)
plot(dnorm(sort(X),0,1),type="l") 
plot(sort(X), pnorm(sort(X), 0,1), type="l", ylab="F(x)", main="Normal(0,1) CDF")

Y<-rnorm(1000,0,1)
summary(Y)
plot(Y)
plot(dnorm(sort(Y),0,1),type="l") 
plot(sort(Y), pnorm(sort(Y), 0,1), type="l", ylab="F(Y)", main="Normal(0,1) CDF")

X<-rnorm(100,2,2)
erro<-2*X+X^2
Y<-X+erro

par(mfrow=c(1,2))
ts.plot(X)
ts.plot(Y)

# Exercício

1) Gere uma amostra da distribuição Poisson de tamanho 100  e média 2 (Use rpois)
2)  Obtenha um resumo descritivo desta amostra
3)  Faça o gráfico da função densidade (Use dpois) e função de distribuição 
(Use ppois) desta amostra.

#Resposta
X<-rpois(100,2)
summary(X)
plot(dpois(sort(X),2)) 
plot(sort(X), ppois(sort(X), 2), type="s", ylab="F(x)", main="Poisson(1) CDF")

#Banco de dados do R
data()
data(mtcars)  
head(mtcars)

dados1 <- matrix(scan('http://verde.esalq.usp.br/~vitor/ST/Dados/Morettin/D-BANESPA.txt'),ncol=1)
dados1<-read.table('C:/Documents and Settings/usuario/Meus documentos/D-BANESPA.txt')
dados_teste2<-scan('C:/Documents and Settings/usuario/Meus documentos/D-BANESPA.txt')
dados2 <- read.csv('http://verde.esalq.usp.br/~vitor/ST/Dados/Morettin/BEBIDA.csv')
dados2
dados1
dados3 <- c(23, 34, 56, 55, 43, 22, 39, 50)
mean(dados1)
median(dados1)
var(dados1)
sd(dados1)
min(dados1)
max(dados1)
sum(dados1)
prod(dados1)

summary(dados1)

data()
data(airquality)
head(airquality)
par(mfrow=c(2,2))
sfit <- lm(airquality$Ozone~airquality$Solar.R)
summary(sfit)
names(sfit)
plot(airquality$Ozone,airquality$Solar.R)
names(sfit)
sfit$coef
abline(sfit$coeff)
title('Dados e Ajuste')
hist(sfit$resid)
title('Residual Distribution')
plot(sfit$fitted,sfit$resid)
title('Residual/Fitted Values')
par(mfrow=c(1,1))

###################################################################
# Sexta aula
###################################################################
# Distribuições Normal

?dnorm
x <- rnorm(1000,2,3);x
hist(x, freq=F)
curve(dnorm(x,2,3), col=2, lty=1, lwd=2, add=T)

?rnorm
# dnorm(q,mean=0,sd=1) Densidade de probabilidade f(x) em x
# pnorm(q,mean=0,sd=1) Função de densidade acumulada F(x) Ex: P(X <= -1)
# qnorm(p,mean=0,sd=1) Quantil para a probabilidade p Ex: P(X <= a) = 0,975
# rnorm(n,mean=0,sd=1) Gera n números - retira uma amostra da distribuição

# f(x) = 1/(???(2??)??) e^-((x - ??)^2/(2??^2))
# quando x = -1, então (1/sqrt(2*pi)) * exp((-1/2)*(-1)^2)
# portanto: 
1/sqrt(2*pi)*exp(-1/2)

# que é o mesmo que:
dnorm(-1)

# gráficos da fdp e fda com + / - 3 desvios padrões da média no intervalo -3 e 3
plot(dnorm, -3, 3)
plot(pnorm, -3, 3)

# alterando as médias e variâncias
plot(function(x) dnorm(x, 100, 8), 60, 140, ylab='f(x)')
plot(function(x) dnorm(x, 90, 8), 60, 140, add=T, col=2)
plot(function(x) dnorm(x, 100, 15), 60, 140, add=T, col=3)
legend(120, 0.05, c("N(100,64)","N(90,64)","N(100,225)"), fill=1:3)

# gerando 100 observações com média 0 e variância 1
X <- rnorm(100,0,1)
summary(X)
plot(X)
plot(dnorm(sort(X),0,1),type="l") 
plot(sort(X), pnorm(sort(X), 0,1), type="l", ylab="F(x)", main="Normal(0,1) FDA")

# gerando 1000 observações com média 0 e variância 1
Y <- rnorm(1000,0,1)
summary(Y)
plot(Y)
plot(dnorm(sort(Y),0,1),type="l") 
plot(sort(Y), pnorm(sort(Y), 0,1), type="l", ylab="F(Y)", main="Normal(0,1) CDF")

# gerando 10000 observações com média 0 e variância 1
Z <- rnorm(10000,0,1)
summary(Z)
plot(Z)
plot(dnorm(sort(Z),0,1),type="l") 
plot(sort(Z), pnorm(sort(Z), 0,1), type="l", ylab="F(Z)", main="Normal(0,1) CDF")

# retormando os exercícios da aula anterior

# P(Z <= 0,32)
pnorm(0.32,0,1)

# P(0 < Z <= 1,71)
pnorm(1.71,0,1)-0.5

# P(1,32 < Z <= 1,79)
pnorm(1.79,0,1)-pnorm(1.32,0,1)

# P(Z >= 1,5)
1 - pnorm(1.5,0,1)

# P(Z <= -1,3)
pnorm(-1.3,0,1)

# P(-1,5 <= Z <= 1,5)
pnorm(1.5,0,1)-pnorm(-1.5,0,1)

# P(-2,3 <= Z <= -1,49)
pnorm(-1.49,0,1)-pnorm(-2.3,0,1)

# P(-1 <= Z <= 2)
pnorm(2,0,1)-pnorm(-1,0,1)

####################################

# P(Z <= z) = 0,975
qnorm(0.975)

# P(0 < Z <= z) = 0,4975
qnorm(0.9975)

# P(Z >= z) = 0.3
qnorm(0.7)

# P(Z >= z) = 0,975
qnorm(1-0.975)

# P(Z <= z) = 0,1
qnorm(0.10)

# P(-z <= Z <= z) = 0,8
qnorm(0.9)

####################################

# P(6 <= X <= 12)
pnorm(12,10,8)-pnorm(6,10,8)

# P(X <= 8 ou X > 14)
pnorm(8,10,8)+(1-pnorm(14,10,8))

####################################

# distribuição t-Student

pt(1.86,8) # = 0,95
qt(0.95,8) # = 1,86

####################################

# criação de uma função
f <- function(x){
  sqrt(x)
}
f(1:4)

gmedia <- function(x){
  prx <- sum(x)
  n <- length(x)
  gm <- prx/n
  gm
}

gmedia(1:10)
mean(1:10)

# crie as seguintes funções e calcule f(x) nos respectivos intervalos. Plote seu gráfico na forma de linha

# a) f(x) = 2x + 3, 0 < x 6
# b) f(x) = 10x - 9, -3 < x 3
# c) f(x) = x^2 + 3, 0 < x 100
# d) f(x) = x^3 + 3x, -20 < x 20
# e) f(x) = x^2 / 2 + 3x - 10, 1000 < x 6000
# f) f(x) = pi/3 + 2*x3 + x^4 / 5


fxa <- function(x){
  2*x + 3
}

ya <- fxa(seq(0,6,0.1))
xa <- seq(0,6,0.1)
plot(xa,ya, type = 'l', col=2, lwd=3)

# Integral de uma dimensão
fx <- function(x) x^2 
integrate(fx, -1, 1)

# importando dados
ICV <- read.table("C:\\Dados\\icv.txt", head=F)
ICV
summary(ICV)
sqrt(var(ICV))
dim(ICV)

is.data.frame(ICV)

icv <- as.vector(ICV[,1])
icv
length(icv)
mean(icv)

calcule as principais estatísticas descritivas
plote seu gráfico
calcule a probabilidade de 80 < ICV < 150
calcule a probabilidade de 100 < ICV < 750
calcule a probabilidade de 230 < media_ICV < 270

# calculando passo a passo a prob. do limite inferior (LI)
media <- mean(icv)
dp <- sqrt(var(icv));dp
dp/sqrt(n)
n <- 114
esquerda <- (230-mean(icv))/(sqrt(var(icv))/(sqrt(n)));esquerda
pnorm(esquerda,0,1)

# montando uma função para calcular a prob do LI
prob_e <- function(x){
  media <- mean(icv)
  dp <- sqrt(var(icv))
  n <- 114
  esquerda <- (x-mean(icv))/(sqrt(var(icv))/(sqrt(n)))
  pnorm(esquerda,0,1)
}
prob_e(230)

# montando uma função para calcular a prob do LS
prob_d <- function(y){
  media <- mean(icv)
  dp <- sqrt(var(icv))
  n <- 114
  direita <- (y+mean(icv))/sqrt(var(icv))/(sqrt(n))
  pnorm(direita,0,1)
}
prob_d(270)

prob_d(270)-prob_e(230)

# montando uma função para calcular a prob final
prob <- function(x,y){
  media <- mean(icv)
  dp <- sqrt(var(icv))
  n <- 114
  esquerda <- (x-mean(icv))/(sqrt(var(icv))/(sqrt(n)))
  prob_es <- pnorm(esquerda,0,1)
  direita <- (y+mean(icv))/sqrt(var(icv))/(sqrt(n))
  prob_di <- pnorm(direita,0,1)
  prob_di-prob_es
}
prob(230,270)


###################################################################
# Sétima aula 
###################################################################
O programa R inclui funcionalidade para operações com distribuições de probabilidades. Para cada distribuição há 4 operações básicas indicadas pelas letras:

d* > calcula a densidade de probabilidade f(x) no ponto
p* > calcula a função de probabilidade acumulada F(x) no ponto
q* > calcula o quantil correspondente a uma dada probabilidade
r* > retira uma amostra da distribuição

Distribuição Uniforme Discreta

Não há entre as funções básicas do R uma função específica para a distribuição uniforme discreta, embora algumas a função sample() seja útil nesse caso:
A função sample() não é restrita à distribuição uniforme discreta, podendo ser usada para sorteios, com ou sem reposição (argumento replace, default sem reposição), com a possibilidade de associar diferentes probabilidades a cada elemento (argumento prob, default probabilidades iguais para os elementos).

args(sample)
function (x, size, replace = FALSE, prob = NULL) 

Vejamos alguns exemplos:

sample(1:5, 10, rep = T)

sorteio de 3 números entre os inteiros de 0 a 20
sample(0:60, 6)

sorteio de 5 números entre os elementos de um certo vetor
x <- c(1, 2, 6, 3, 10) 
sample(x, 2)

sorteio de 100 números entre os possíveis resultados do lançamento de um dado, com reposição
sample(1:6, 100, rep = T)


Distribuição Binomial

Seja X uma v.a. com distribuição Binomial com n = 10 e p = 0.35. Vamos ver os comandos do R para:

fazer o gráfico das função de densidade
idem para a função de probabilidade
calcular P [X = 7]
calcular P [X < 8] = P[X = 7]
calcular P [X = 8] = P[X > 7]
calcular P [3 < X = 6] = P[4 = X < 7]

x <- 0:10 
fx <- dbinom(x, 10, 0.35) 
plot(x, fx, type = "h") 
Fx <- pbinom(x, 10, 0.35) 
plot(x, Fx, type = "s")
dbinom(7, 10, 0.35)
pbinom(7, 10, 0.35)
sum(dbinom(0:7, 10, 0.35))
1 - pbinom(7, 10, 0.35)
pbinom(6, 10, 0.35) - pbinom(3, 10, 0.35)
sum(dbinom(4:6, 10, 0.35))

Distruição Poisson

dpois(x, lambda, log = FALSE)
ppois(q, lambda, lower.tail = TRUE, log.p = FALSE)
qpois(p, lambda, lower.tail = TRUE, log.p = FALSE)
rpois(n, lambda)

Arguments

x - vector of (non-negative integer) quantiles.
q - vector of quantiles.
p - vector of probabilities.
n - number of random values to return.
lambda - vector of (non-negative) means.
lower.tail - logical; if TRUE (default), probabilities are P[X = x], otherwise, P[X > x].


# Exercício

1) Gere uma amostra da distribuição Poisson de tamanho 100  e média 2 (Use rpois)
2)  Obtenha um resumo descritivo desta amostra
3)  Faça o gráfico da função densidade (Use dpois) e função de distribuição 
(Use ppois) desta amostra.

#Resposta
X<-rpois(100,2)
summary(X)
plot(dpois(sort(X),2)) 
plot(sort(X), ppois(sort(X), 2), type="s", ylab="F(x)", main="Poisson(1) CDF")