In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Descricao
#
# Regiao de confianca (elipse) - dados sobre radiacao em microondas
#
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(tseries)) # Para o teste de normalidade de Jarque-Bera
silence(library(ellipse))
In [2]:
#==============================================================================
# Lendo os dados
#==============================================================================
x <- read.table("../dados/d03_microondas.txt",header=T)
print(head(x))
Fechada Aberta 1 0,15 0,30 2 0,09 0,09 3 0,18 0,30 4 0,10 0,10 5 0,05 0,10 6 0,12 0,12
In [3]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
hist(x[,1],prob=T,col="red",main="",xlab="Fechada",ylab="Densidade")
hist(x[,2],prob=T,col="red",main="",xlab="Aberta",ylab="Densidade")
In [4]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
boxplot(x[,1],col="red",main="",xlab="Fechada")
boxplot(x[,2],col="red",main="",xlab="Aberta")
In [5]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
qqnorm(x[,1],col="red",main="Fechada",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3))
qqline(x[,1],lwd=2,col="darkgreen")
qqnorm(x[,2],col="red",main="Aberta",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3))
qqline(x[,2],lwd=2,col="darkgreen")
In [6]:
#==============================================================================
# Testando normalidade: teste de Jarque-Bera (parametrico)
#==============================================================================
jarque.bera.test(x[,1])
jarque.bera.test(x[,2])
Jarque Bera Test data: x[, 1] X-squared = 11,939, df = 2, p-value = 0,002555
Jarque Bera Test data: x[, 2] X-squared = 20,969, df = 2, p-value = 2,797e-05
In [7]:
#==============================================================================
# Testando normalidade: teste de Shapiro-Wilk (nao parametrico)
#==============================================================================
shapiro.test(x[,1])
shapiro.test(x[,2])
Shapiro-Wilk normality test data: x[, 1] W = 0,85793, p-value = 9,902e-05
Shapiro-Wilk normality test data: x[, 2] W = 0,8292, p-value = 1,976e-05
In [8]:
#==============================================================================
# Transformacao dos dados para tentar se adequar a hipotese de normalidade
#==============================================================================
z <- x^0.25
In [9]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
hist(z[,1],prob=T,col="darkorange",main="",xlab="Fechada",ylab="Densidade")
hist(z[,2],prob=T,col="darkorange",main="",xlab="Aberta",ylab="Densidade")
In [10]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
boxplot(z[,1],col="darkorange",main="",xlab="Fechada")
boxplot(z[,2],col="darkorange",main="",xlab="Aberta")
In [11]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(1,2),bty="n")
qqnorm(z[,1],col="darkorange",main="Fechada",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3))
qqline(z[,1],lwd=2,col="darkgreen")
qqnorm(z[,2],col="darkorange",main="Aberta",pch=15,xlab="Percentil teórico",
ylab="Percentil amostral",xlim=c(-3,3))
qqline(z[,2],lwd=2,col="darkgreen")
In [12]:
#==============================================================================
# Testando normalidade: teste de Jarque-Bera (parametrico)
#==============================================================================
library(tseries)
jarque.bera.test(z[,1])
jarque.bera.test(z[,2])
Jarque Bera Test data: z[, 1] X-squared = 0,21782, df = 2, p-value = 0,8968
Jarque Bera Test data: z[, 2] X-squared = 0,10021, df = 2, p-value = 0,9511
In [13]:
#==============================================================================
# Testando normalidade: teste de Shapiro-Wilk (nao parametrico)
#==============================================================================
shapiro.test(z[,1])
shapiro.test(z[,2])
Shapiro-Wilk normality test data: z[, 1] W = 0,96481, p-value = 0,2188
Shapiro-Wilk normality test data: z[, 2] W = 0,94282, p-value = 0,03589
In [14]:
#==============================================================================
# Regiao (elipse) de confianca.
#==============================================================================
alpha <- 0.05
n <- dim(z)[1]
p <- dim(z)[2]
xbarra <- apply(z,2,"mean")
S <- var(z)
IS <- solve(S)
DES <- eigen(S)
c2 <- (n-1)*p*qf(1-alpha,p,n-p)/(n-p)
eval <- DES$val
evec <- DES$vec
semiz1 <- sqrt(c2*eval[1]/n)
semiz2 <- sqrt(c2*eval[2]/n)
y1 <- xbarra + semiz1*evec[,1] # definindo semi-eixo 1
y2 <- xbarra + semiz2*evec[,2] # definindo semi-eixo 2
y3 <- xbarra - semiz1*evec[,1] # definindo semi-eixo 1
y4 <- xbarra - semiz2*evec[,2] # definindo semi-eixo 2
In [15]:
#==============================================================================
# Codigo para grafico da elipse de confianca
# Note que o tamanho da amostra n esta' divindo c2, isto e', c2/n
#==============================================================================
v <- ellipse(S,t=sqrt(c2/n),centre=xbarra,npoints=2000)
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(4.5,4.5,2,2),xpd=T,bty="n")
plot(v,col="darkgreen",cex=1.2,type="l",lwd=3,xlab="Porta fechada",
ylab="Porta aberta",xlim=c(0.51,0.63),ylim=c(0.54,0.66))
arrows(xbarra[1],xbarra[2],y1[1],y1[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y2[1],y2[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y3[1],y3[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y4[1],y4[2],col="darkblue",lwd=2,code=0,lty=2)
In [16]:
#==============================================================================
# Intervalos T^2
#==============================================================================
lim.inf <- xbarra - sqrt(c2)*sqrt(diag(S)/n)
lim.sup <- xbarra + sqrt(c2)*sqrt(diag(S)/n)
lim <- cbind(lim.inf,lim.sup)
print(lim.inf)
print(lim.sup)
Fechada Aberta 0,5166803 0,5550817 Fechada Aberta 0,6118347 0,6508807
In [17]:
#==============================================================================
# Elipse de confianca e intervalo T^2
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(4.5,4.5,2,2),bty="n")
plot(v,col="darkgreen",cex=1.2,type="l",lwd=3,xlab="Porta fechada",
ylab="Porta aberta",xlim=c(0.51,0.63),ylim=c(0.54,0.66))
arrows(xbarra[1],xbarra[2],y1[1],y1[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y2[1],y2[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y3[1],y3[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y4[1],y4[2],col="darkblue",lwd=2,code=0,lty=2)
abline(v=lim[1,],lwd=1,col="black",lty=3)
abline(h=lim[2,],lwd=1,col="black",lty=3)
lines(rep(lim[1,1],2),lim[2,],lwd=3,col="red",lty=2)
lines(rep(lim[1,2],2),lim[2,],lwd=3,col="red",lty=2)
lines(lim[1,],rep(lim[2,1],2),lwd=3,col="red",lty=2)
lines(lim[1,],rep(lim[2,2],2),lwd=3,col="red",lty=2)
In [18]:
#==============================================================================
# Intervalos via Bonferroni
#==============================================================================
tlim.inf <- xbarra - qt(1-alpha/(2*p),n-1)*sqrt(diag(S)/n)
tlim.sup <- xbarra + qt(1-alpha/(2*p),n-1)*sqrt(diag(S)/n)
tlim <- cbind(tlim.inf,tlim.sup)
print(tlim.inf)
print(tlim.sup)
Fechada Aberta 0,5212495 0,5596819 Fechada Aberta 0,6072655 0,6462806
In [19]:
#==============================================================================
# Elipse de confianca, intervalo T^2 e intervalo de Bonferroni
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,5),
mar=c(4.5,4.5,2,2),bty="n")
plot(v,col="darkgreen",cex=1.2,type="l",lwd=3,xlab="Porta fechada",
ylab="Porta aberta",xlim=c(0.51,0.63),ylim=c(0.54,0.66))
arrows(xbarra[1],xbarra[2],y1[1],y1[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y2[1],y2[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y3[1],y3[2],col="darkblue",lwd=2,code=0,lty=2)
arrows(xbarra[1],xbarra[2],y4[1],y4[2],col="darkblue",lwd=2,code=0,lty=2)
abline(v=lim[1,],lwd=1,col="black",lty=3)
abline(h=lim[2,],lwd=1,col="black",lty=3)
lines(rep(lim[1,1],2),lim[2,],lwd=3,col="red",lty=2)
lines(rep(lim[1,2],2),lim[2,],lwd=3,col="red",lty=2)
lines(lim[1,],rep(lim[2,1],2),lwd=3,col="red",lty=2)
lines(lim[1,],rep(lim[2,2],2),lwd=3,col="red",lty=2)
abline(v=tlim[1,],lwd=1,col="black",lty=3)
abline(h=tlim[2,],lwd=1,col="black",lty=3)
lines(rep(tlim[1,1],2),tlim[2,],lwd=3,col="magenta",lty=2)
lines(rep(tlim[1,2],2),tlim[2,],lwd=3,col="magenta",lty=2)
lines(tlim[1,],rep(tlim[2,1],2),lwd=3,col="magenta",lty=2)
lines(tlim[1,],rep(tlim[2,2],2),lwd=3,col="magenta",lty=2)
