quinta-feira, 15 de outubro de 2009

Estrutura do arquivo cereal3.txt

Códigos em R para exercício 4 , lista 3 - A.M.

## observacao:
## o arquivo cereal3.txt deve estar salvo
## no mesmo diretorio onde o R está,
## ou mudar o diretorio do R para o diretorio
## onde está salvo o arquivo cereal3.txt


x<- read.table("cereal3.txt", head=T) # lendo os dados
x1<- x[,1]
x1<- t(t(x1))
x2<- x[,2]
x2<- t(t(x2))
x3<- x[,3]
x3<- t(t(x3))
x4<- x[,4]
x4<- t(t(x4))
X <- cbind(x1,x2,x3,x4) # matriz de dados

mean.x1<- mean(x1) # media e dp. das variaveis
mean.x2<- mean(x2)
mean.x3<- mean(x3)
mean.x4<- mean(x4)
dp.x1<- sqrt(var(x1))
dp.x2<- sqrt(var(x2))
dp.x3<- sqrt(var(x3))
dp.x4<- sqrt(var(x4))
z1<- (X[,1]- mean.x1)/ dp.x1
z1<- t(t(z1))
z2<- (X[,2]- mean.x2)/ dp.x2
z2<- t(t(z2))
z3<- (X[,3]- mean.x3)/ dp.x3
z3<- t(t(z3))
z4<- (X[,4]- mean.x4)/ dp.x4
Z<- cbind(z1,z2,z3,z4) # matriz de dados padronizada

m.cor<-cor(Z)
p <- ncol(Z)
aut.val <- eigen(m.cor)$values
aut.vec <- -(eigen(m.cor)$vectors)
m.aut.val <- t(matrix((sqrt(aut.val)),p,p))
result.cp.cor <- princomp(Z,cor=TRUE)
corr.cp.var <- aut.vec*m.aut.val
screeplot(result.cp.cor,type=c("lines"),main="autovalores") #grafico screeplot
comp.princ <- Z %*% aut.vec
cp1<- comp.princ[,1]
cp1<- t(t(cp1))
cp2<- comp.princ[,2]
cp2<- t(t(cp2))
plot(cp1,cp2)
qqnorm(cp1)
qqline(cp1)

Códigos em R para exercício 5 , lista 3 - A.M.

x <- iris # matriz de dados completa
x1<- x[,1]
x1<- t(t(x1))
x2<- x[,2]
x2<- t(t(x2))
x3<- x[,3]
x3<- t(t(x3))
x4<- x[,4]
x4<- t(t(x4))
X <- cbind(x1,x2,x3,x4) # matriz de dados
mean.x1<- mean(x1) # media e dp. das variaveis
mean.x2<- mean(x2)
mean.x3<- mean(x3)
mean.x4<- mean(x4)
dp.x1<- sqrt(var(x1))
dp.x2<- sqrt(var(x2))
dp.x3<- sqrt(var(x3))
dp.x4<- sqrt(var(x4))
z1<- (X[,1]- mean.x1)/ dp.x1
z1<- t(t(z1))
z2<- (X[,2]- mean.x2)/ dp.x2
z2<- t(t(z2))
z3<- (X[,3]- mean.x3)/ dp.x3
z3<- t(t(z3))
z4<- (X[,4]- mean.x4)/ dp.x4
Z<- cbind(z1,z2,z3,z4) # matriz de dados padronizada
nomes<-cbind("Sepal.Length"," Sepal.Width", "Petal.Length"," Petal.Width")
colnames(Z)<- nomes # nomeando as colunas
m.X.completa<- Z #m.X.completa recebe a matriz de dados padronizada
m.X.completa<- Z
v.grupos<- matrix(nrow=150,ncol=1)
v.grupos[1:50,1]<-1
v.grupos[51:100,1]<-2
v.grupos[101:150,1]<-3
G<- 3
v.n<- matrix(nrow=3,ncol=1)
v.n[1,1]<-50
v.n[2,1]<-50
v.n[3,1]<-50



### aproveitando o codigo do prof, para calcular a matriz Sp^2
grupo <- 1
p<- ncol(m.X.completa)
m.X.k <- m.X.completa[v.grupos==grupo,]
Sigma.k <- cov(m.X.k)
m.Sigma.completa <- cbind(grupo,Sigma.k)
Sigma.P <- (v.n[grupo]-1)*Sigma.k # estimativa ponderada
aux.k.1 <- (v.n[grupo] - 1)*log(det(Sigma.k))
grupo <- grupo + 1
for (i in 2:G)
{
m.X.k <- m.X.completa[v.grupos==grupo,] # pegar os dados referentes ao grupo i
Sigma.k <- cov(m.X.k)
m.Sigma.completa <- rbind(m.Sigma.completa,cbind(grupo,Sigma.k))
Sigma.P <- Sigma.P + (v.n[grupo]-1)*Sigma.k # estimativa ponderada
aux.k.1 <- aux.k.1 + (v.n[grupo] - 1)*log(det(Sigma.k))
grupo <- grupo + 1
}

sigma1<- sqrt(Sigma.P[1,1])
sigma2<- sqrt(Sigma.P[2,2])
sigma3<- sqrt(Sigma.P[3,3])
sigma4<- sqrt(Sigma.P[4,4])

D<-matrix(nrow=4,ncol=4) # so gambiarra, to manjando!
D[1,1]<- sigma1
D[2,2]<- sigma2
D[3,3]<- sigma3
D[4,4]<- sigma4
D<- diag(D)
D<- diag(D) # D é matriz diagonal de desvios padroes
R<- solve(D) %*% Sigma.P %*% solve(D) # matriz de corr. amostral

aut.val <- eigen(R)$values
aut.vec <- -(eigen(R)$vectors)
m.aut.val <- t(matrix((sqrt(aut.val)),p,p))
result.cp.cor <- princomp(Z,cor=TRUE)
corr.cp.var <- aut.vec*m.aut.val
summary(result.cp.cor)
screeplot(result.cp.cor,type=c("lines"),main="autovalores")

comp.princ <- Z %*% aut.vec # componentes principais