In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library(tseries)) # Para o teste de normalidade de Jarque-Bera
silence(library(MASS))
In [2]:
#==============================================================================
# Lendo dados em d08_lagarto.txt
#==============================================================================
x <- matrix(scan("../dados/d08_lagarto.txt"),,3,T)
gen <- c(0,1,0,0,1,0,1,0,1,0,1,0,1,1,1,1,0,1,1,1,0,0,1,0,0)
# f m f f m f m f m f m f m m m m f m m m f f m f f
massa <- x[,1]
svl <- x[,2]
hls <- x[,3]
In [3]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(2,2),bty="n",lab=c(6,6,5))
hist(massa,prob=T,col="red",main="",xlab="Massa",ylab="Densidade",
xlim=c(0,16),ylim=c(0,0.2))
hist(svl, prob=T,col="red",main="",xlab="SVL",ylab="Densidade",
xlim=c(40,100),ylim=c(0,0.06))
hist(hls, prob=T,col="red",main="",xlab="HLS",ylab="Densidade",
xlim=c(80,180),ylim=c(0,0.035))
In [4]:
#==============================================================================
# Verificando normalidade
#==============================================================================
par(mfrow=c(2,2),lab=c(6,6,5),bty="n")
boxplot(massa,col="red",main="",xlab="Massa",ylim=c(0,16))
boxplot(svl, col="red",main="",xlab="SVL",ylim=c(40,100))
boxplot(hls, col="red",main="",xlab="HLS",ylim=c(80,180))
In [5]:
#==============================================================================
# Testando normalidade
# Teste de Shapiro-Wilk
#==============================================================================
shapiro.test(massa)
shapiro.test(svl)
shapiro.test(hls)
Shapiro-Wilk normality test data: massa W = 0,97981, p-value = 0,8814
Shapiro-Wilk normality test data: svl W = 0,97248, p-value = 0,7083
Shapiro-Wilk normality test data: hls W = 0,97147, p-value = 0,6824
In [6]:
#==============================================================================
# Testando normalidade
# Teste de Jarque-Bera
#==============================================================================
jarque.bera.test(massa)
jarque.bera.test(svl)
jarque.bera.test(hls)
Jarque Bera Test data: massa X-squared = 0,55969, df = 2, p-value = 0,7559
Jarque Bera Test data: svl X-squared = 1,1667, df = 2, p-value = 0,558
Jarque Bera Test data: hls X-squared = 0,15422, df = 2, p-value = 0,9258
In [7]:
#==============================================================================
# Analise de discriminante - probabilidades a priori proporcionais
#==============================================================================
disc1 <- lda( gen ~ massa + svl + hls )
print(disc1)
Call:
lda(gen ~ massa + svl + hls)
Prior probabilities of groups:
0 1
0,48 0,52
Group means:
massa svl hls
0 7,011583 63,04167 118,0000
1 10,232769 73,34615 139,7692
Coefficients of linear discriminants:
LD1
massa -0,57233549
svl -0,09082882
hls 0,29489697
In [8]:
#==============================================================================
# Erros de classificacao
#==============================================================================
z <- c(predict(disc1,data.frame(cbind(massa,svl,hls)))$class)
print(as.numeric(z)-1)
# f m f f m f m f m f m f m m m m f m m m f f m f f
[1] 0 1 0 0 1 1 1 0 1 0 1 0 1 1 1 1 0 1 1 1 0 0 1 0 0
In [9]:
#==============================================================================
# Machos
#==============================================================================
print(sum(gen)) # machos
print(sum(as.numeric(z)-1)) # classificados como machos
[1] 13 [1] 14
In [10]:
#==============================================================================
# Femeas
#==============================================================================
print(sum(1-gen)) # femea
print(sum(2-as.numeric(z))) # classificados como femea
[1] 12 [1] 11
In [11]:
#==============================================================================
# Analise de discriminante - probabilidades a priori iguais
#==============================================================================
disc2 <- lda( gen ~ massa + svl + hls, prior=c(0.5,0.5) )
print(disc2)
Call:
lda(gen ~ massa + svl + hls, prior = c(0.5, 0.5))
Prior probabilities of groups:
0 1
0,5 0,5
Group means:
massa svl hls
0 7,011583 63,04167 118,0000
1 10,232769 73,34615 139,7692
Coefficients of linear discriminants:
LD1
massa -0,57233549
svl -0,09082882
hls 0,29489697
In [12]:
#==============================================================================
# Erros de classificacao
#==============================================================================
w <- c(predict(disc2,data.frame(cbind(massa,svl,hls)))$class)
print(as.numeric(w)-1)
# f m f f m f m f m f m f m m m m f m m m f f m f f
[1] 0 1 0 0 1 0 1 0 1 0 1 0 1 1 1 1 0 1 1 1 0 0 1 0 0
In [13]:
#==============================================================================
# Machos
#==============================================================================
print(sum(gen)) # machos
print(sum(as.numeric(w)-1)) # classificados como machos
[1] 13 [1] 13
In [14]:
#==============================================================================
# Femeas
#==============================================================================
print(sum(1-gen)) # femea
print(sum(2-as.numeric(w))) # classificados como femea
[1] 12 [1] 12