In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo para gerar de uma normal multivariada
#==============================================================================
#
# Parte 1: utilizando a decomposicao de Cholesky
#
#==============================================================================
set.seed(54345) # Semente
n <- 100000 # Tamanho da amostra
#==============================================================================
MU <- c(-3,3,1)
SIGMA <- matrix(c(1,0.9,0.8,0.9,1,0.7,0.8,0.7,1),3,3,T)
In [2]:
#==============================================================================
# Inicio (Cholesky)
#==============================================================================
U <- chol(SIGMA)
print("===============================")
print(t(U))
print("===============================")
print(t(U)%*%U) # Conferindo SIGMA
print("===============================")
print(SIGMA)
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,00000000 0,000000
[2,] 0,9 0,43588989 0,000000
[3,] 0,8 -0,04588315 0,598243
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
In [3]:
#==============================================================================
# Gerando da normal multivariada (trivariada) via Cholesky
#==============================================================================
y <- matrix(NA,n,3)
for(i in 1:n)
y[i,] <- MU + t(U)%*%rnorm(3)
print("===============================")
print(apply(y,2,"mean")) #(-3,3,1)
print("===============================")
print(cov(y))
print("===============================")
print(cor(y))
print("===============================")
print(SIGMA)
[1] "==============================="
[1] -2,996126 3,003515 1,002702
[1] "==============================="
[,1] [,2] [,3]
[1,] 0,9981761 0,9012971 0,7995698
[2,] 0,9012971 1,0027973 0,7004597
[3,] 0,7995698 0,7004597 0,9990415
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0000000 0,9008611 0,8006838
[2,] 0,9008611 1,0000000 0,6998175
[3,] 0,8006838 0,6998175 1,0000000
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
In [4]:
#==============================================================================
# Dispersão
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(5,5,5),
mar=c(4.5,5,2,1),bty="n")
colnames(y) <- c(expression(y[1]),expression(y[2]),expression(y[3]))
pairs(y,pch=".")
In [5]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[1]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,1],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[1]),
ylab=expression(f(y[1])))
sig1 <- sqrt(SIGMA[1,1])
xseq <- seq(MU[1]-4*sig1,MU[1]+4*sig1,length=1000)
yseq <- dnorm(xseq,MU[1],sig1)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,1],main="",col="darkgreen")
qqnorm(y[,1],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,1],lwd=3)
In [6]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[2]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,2],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[2]),
ylab=expression(f(y[2])))
sig2 <- sqrt(SIGMA[2,2])
xseq <- seq(MU[2]-4*sig2,MU[2]+4*sig2,length=1000)
yseq <- dnorm(xseq,MU[2],sig2)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,2],main="",col="darkgreen")
qqnorm(y[,2],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,2],lwd=3)
In [7]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[3]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,3],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[3]),
ylab=expression(f(y[3])))
sig2 <- sqrt(SIGMA[3,3])
xseq <- seq(MU[3]-4*sig2,MU[3]+4*sig2,length=1000)
yseq <- dnorm(xseq,MU[3],sig2)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,2],main="",col="darkgreen")
qqnorm(y[,2],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,2],lwd=3)
In [8]:
#==============================================================================
#
# Parte 2: utilizando a "SVD"
#
#==============================================================================
rm(list=ls())
set.seed(54345) # Semente
n <- 100000 # Tamanho da amostra
#==============================================================================
MU <- c(-3,3,1)
SIGMA <- matrix(c(1,0.9,0.8,0.9,1,0.7,0.8,0.7,1),3,3,T)
In [9]:
#==============================================================================
# Inicio (SVD)
#==============================================================================
decomp <- svd(SIGMA)
V <- decomp$u # autovetores
lambda <- decomp$d # autovalores
H <- V%*%diag(sqrt(lambda))
print("===============================")
print(H)
print("===============================")
print(H%*%t(H)) # Conferindo SIGMA
print("===============================")
print(SIGMA)
[1] "==============================="
[,1] [,2] [,3]
[1,] -0,9686611 0,1074754 0,22393003
[2,] -0,9330609 0,3173500 -0,16937043
[3,] -0,8909961 -0,4491760 -0,06608266
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
In [10]:
#==============================================================================
# Gerando da normal multivariada (trivariada) via SVD
#==============================================================================
y <- matrix(NA,n,3)
for(i in 1:n)
y[i,] <- MU + H%*%rnorm(3)
print("===============================")
print(apply(y,2,"mean")) #(-3,3,1)
print("===============================")
print(cov(y))
print("===============================")
print(cor(y))
print("===============================")
print(SIGMA)
[1] "==============================="
[1] -3,0038932 2,9965174 0,9965625
[1] "==============================="
[,1] [,2] [,3]
[1,] 0,9953080 0,8951809 0,8014065
[2,] 0,8951809 0,9951111 0,7002443
[3,] 0,8014065 0,7002443 1,0027457
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0000000 0,8994899 0,8021927
[2,] 0,8994899 1,0000000 0,7010006
[3,] 0,8021927 0,7010006 1,0000000
[1] "==============================="
[,1] [,2] [,3]
[1,] 1,0 0,9 0,8
[2,] 0,9 1,0 0,7
[3,] 0,8 0,7 1,0
In [11]:
#==============================================================================
# Dispersão
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(5,5,5),
mar=c(4.5,5,2,1),bty="n")
colnames(y) <- c(expression(y[1]),expression(y[2]),expression(y[3]))
pairs(y,pch=".")
In [12]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[1]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,1],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[1]),
ylab=expression(f(y[1])))
sig1 <- sqrt(SIGMA[1,1])
xseq <- seq(MU[1]-4*sig1,MU[1]+4*sig1,length=1000)
yseq <- dnorm(xseq,MU[1],sig1)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,1],main="",col="darkgreen")
qqnorm(y[,1],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,1],lwd=3)
In [13]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[2]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,2],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[2]),
ylab=expression(f(y[2])))
sig2 <- sqrt(SIGMA[2,2])
xseq <- seq(MU[2]-4*sig2,MU[2]+4*sig2,length=1000)
yseq <- dnorm(xseq,MU[2],sig2)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,2],main="",col="darkgreen")
qqnorm(y[,2],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,2],lwd=3)
In [14]:
#==============================================================================
# Histograma, boxplot e qqnorm para y[3]
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
hist(y[,3],nclass=50,prob=T,main="",col="darkgreen",xlab=expression(y[3]),
ylab=expression(f(y[3])))
sig2 <- sqrt(SIGMA[3,3])
xseq <- seq(MU[3]-4*sig2,MU[3]+4*sig2,length=1000)
yseq <- dnorm(xseq,MU[3],sig2)
lines(xseq,yseq,col="red",lwd=3)
boxplot(y[,2],main="",col="darkgreen")
qqnorm(y[,2],main="",col="darkgreen",
xlab="Percentis teóricos",ylab="Percentis amostrais",pch=16)
qqline(y[,2],lwd=3)
In [15]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================