CHP7266 - Estatástica e Experimentação (em R)

---

###Aula 5 Prof. Lorenzo Zanette - l.zanette@ufc.br

---

---

Distribuições

Alguns exemplos no R:

par(mfrow=c(2,3))
hist(rnorm(1000,5,1),col="red", main="Normal")
hist(rlnorm(1000,5,1),col="orange",main="logNormal")
hist(rt(1000,10),col="blue",main="T de GUINNESS")
hist(rgamma(1000,5),col="red",main="Gamma")
hist(rf(1000,10,10),col="orange",main="F")
hist(rchisq(1000,10),col="blue",main="Chi^2")

---

+ Explorações gráficas

x<-seq(0,10,0.1)
x

##   [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0  1.1  1.2  1.3  1.4
##  [16]  1.5  1.6  1.7  1.8  1.9  2.0  2.1  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9
##  [31]  3.0  3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9  4.0  4.1  4.2  4.3  4.4
##  [46]  4.5  4.6  4.7  4.8  4.9  5.0  5.1  5.2  5.3  5.4  5.5  5.6  5.7  5.8  5.9
##  [61]  6.0  6.1  6.2  6.3  6.4  6.5  6.6  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4
##  [76]  7.5  7.6  7.7  7.8  7.9  8.0  8.1  8.2  8.3  8.4  8.5  8.6  8.7  8.8  8.9
##  [91]  9.0  9.1  9.2  9.3  9.4  9.5  9.6  9.7  9.8  9.9 10.0

par(mfrow=c(1,2))
y<-exp(x) 
plot(y~x,type="l",main="Exponencial")
y<-log(x) 
plot(y~x,type="l",main="Logaritmico")

x<-rnorm(1000,5,2)
y<-exp(x) 
plot(y~x,main="Expo")
y<-log(x) 

## Warning in log(x): NaNs produced

plot(y~x,main="Log")

Quantas variáveis colocar no gráfico?

read.table("legu.txt",header=T)->legu
names(legu)

## [1] "fungi"  "bacter" "pH"

plot(legu$bacter,legu$pH) 

# adiciona texto ao gráfico usando valores em “fungi”
plot(legu$bacter,legu$pH,type="n")
text(legu$bacter,legu$pH,labels=round(legu$fungi,2),pos=1,offset=0.5,cex=0.7)

# alterar cores conforme valores em “FR”
plot(legu$bacter,legu$pH,pch=16,col=ifelse(legu$fungi>2,"red","black"), main="Viva o rubro negro!")

read.table("minhocas.txt",header=T)->minho

names(minho)

## [1] "local"        "area"         "declive"      "vegetacao"    "phsolo"      
## [6] "encharcado"   "densiminhoca"

#checando se a var. é considerada categórica
is.factor(minho$vegetacao)

## [1] FALSE

as.factor(minho$vegetacao)->minho$vegetacao
plot(minho$densiminhoca~minho$vegetacao)

Box plot – resume muita informação – mediana , 25% -75%, variabilidade e possíveis “outliers”

read.table("ozonio.txt",header=T)->ozonio
pairs(ozonio)

coplot(ozonio$ozo~ozonio$vento|ozonio$temp,panel = panel.smooth)

#install.packages("maps")
library(maps)
coplot(lat~long|depth,data= quakes, number=4,panel=function(x,y,...){
usr<-par("usr") 
rect(usr[1],usr[3],usr[2],usr[4],col="light blue") 
map("world2",regions=c("New Zealand", "Fiji"), add=TRUE,lwd=0.1,fill=TRUE,col="dark green")
text(180,-13,"Fiji",adj=1,cex=0.7) 
text(170,-35,"NZ",cex=0.7)
points(x, y,pch="+",cex=0.5)
})

---

---

Buscando a normalidade

***

Teorema do Limite Central : Se amostras repetidas forem retiradas de uma população com variância finita e sua média for calculada, estas médias terão um distribuição normal`

hist(runif(10000,min=1,max=10),main="Distribuição Uniforme",col=terrain.colors(10))

> Ou seja, qualquer valor tem a mesma chance de acontecer…

Se amostrarmos e obtivermos médias das amostras…

means<-numeric(10000)
for (i in 1:10000){ means[i]<-mean(runif(5,min=1,max=10)) }
hist(means,col=terrain.colors(10))

Base para função densidade da Distribuição Normal Padrão

Eq.1 y = exp(-|x|^m)

par(mfrow=c(2,2))
x <- seq(-3,3,0.01)
y <- exp(-abs(x))
plot(x,y,type="l",main= "x")
y <- exp(-abs(x)^2)
plot(x,y,type="l",main= "x^2")
y <- exp(-abs(x)^3)
plot(x,y,type="l",main= "x^3")
y <- exp(-abs(x)^8)
plot(x,y,type="l",main= "x^8")

---

---

---

---

Foram coletados 100 peixes, a média de comprimento = 85cm, d.p = 4cm

Qual a probabilidade de um indivíduo escolhido aleatoriamente seja:

1. Menor do que um comprimento escolhido?

2. Maior do que um comprimento escolhido?

3. Entre dois comprimentos específicos?

4. Menor que 80cm?

Precisamos da Distribuição Normal Padrão, uma DN com média = 0 e d.p = 1

means<-rnorm(100000,0,1)
hist(means,col="light pink",main="DNP")

Quantos d.p. 80cm representa em uma DNP

Ou seja, quanto é 80 cm em zs?

z= (valor - média)/d.p

z<-(80-85)/4

z

## [1] -1.25

# Qual a probabilidade de um valor de -1.25 ou menor em uma DNP?

pnorm(-1.25)

## [1] 0.1056498

Resposta 1 : ca. 10%

---

2. Maior que 92.5cm?

z2<-(92.5-85)/4
z2

## [1] 1.875

# Qual a probabilidade de um valor de 1.875 ou menor em uma DNP?
pnorm(z2)

## [1] 0.9696036

No entanto, pnorm fornece a probabilidade de um valor igual ou menor do que x acontecer…portanto, precisamos de :

1-pnorm(z2)

## [1] 0.03039636

Resposta 2: ca. 3%

---

3. Peixe entre 82.5 e 90 cm?

z3<-(82.5-85)/4
z3

## [1] -0.625

z4<-(90-85)/4
z4

## [1] 1.25

#subtraindo a menor prob. da maior...
pnorm(z4)-pnorm(z3)

## [1] 0.6283647

Resposta 3: 63%

Área da função densidade de probabiliade = 1 = 100%

---

---

Aula 6

Testando uma amostra

alturas <- women[,1]*2.54

hist(alturas,breaks=4,main="",col="light blue")

alturas

##  [1] 147.32 149.86 152.40 154.94 157.48 160.02 162.56 165.10 167.64 170.18
## [11] 172.72 175.26 177.80 180.34 182.88

summary(alturas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   147.3   156.2   165.1   165.1   174.0   182.9

Mas, e a Ana é suas eficientes colegas??? (média = 1.85 m)

qqnorm(alturas)
qqline(alturas)

# Calculando algumas estatísticas conhecidas
desvmed<-mean(alturas)-185
errop<-sd(alturas)/sqrt(15)

#calculando o "t"
tvalor<-desvmed/errop

# qual a prob. de obtermos esse valor com uma amostra de 15 
2*pt(tvalor,14)

## [1] 8.807405e-06

t.test(alturas,mu=185)

## 
##  One Sample t-test
## 
## data:  alturas
## t = -6.785, df = 14, p-value = 8.807e-06
## alternative hypothesis: true mean is not equal to 185
## 95 percent confidence interval:
##  158.8095 171.3905
## sample estimates:
## mean of x 
##     165.1

Outro método: simulando amostras - Bootstrap

resam<-sample(alturas,replace=T)
alturas

##  [1] 147.32 149.86 152.40 154.94 157.48 160.02 162.56 165.10 167.64 170.18
## [11] 172.72 175.26 177.80 180.34 182.88

resam

##  [1] 182.88 170.18 154.94 162.56 180.34 165.10 154.94 177.80 180.34 165.10
## [11] 149.86 152.40 152.40 157.48 147.32

# automatizando 10000 re-amostragens
simu <- numeric(10000)
for(i in 1:10000) simu[i] <- mean(sample(alturas,replace=T))
hist(simu,main="",col="light blue")

summary(simu)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   154.9   163.2   165.1   165.1   167.0   175.6

# em "simu" quantos valores são iguais ou maiores que 185
simu[simu>=185]

## numeric(0)

Ou seja P de que a média seja 185 na pop. em amostras é < 0.000001