In [1]:
# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Exemplo 
#==============================================================================
# A matriz RHO e' positiva definida
#==============================================================================
RHO    <- matrix(c(1,0.9,0.7,0.9,1,0.4,0.7,0.4,1),3,3,T)
print(RHO)

     [,1] [,2] [,3]
[1,]  1,0  0,9  0,7
[2,]  0,9  1,0  0,4
[3,]  0,7  0,4  1,0
In [2]:
#==============================================================================
# Note que e' possivel obter uma solucao via CP, mas ela nao e' exata. 
# So' seria exata para m=p=3
#==============================================================================
dec    <- eigen(RHO)
delta  <- dec$values
P      <- dec$vectors
ell.1  <- sqrt(delta[1])*P[,1]
ell.2  <- sqrt(delta[2])*P[,2]
ell.3  <- sqrt(delta[3])*P[,3]
L      <- cbind(ell.1,ell.2,ell.3)
print(RHO-ell.1%*%t(ell.1))  # Se considerar m=1 a solucao nao e' exata

            [,1]        [,2]        [,3]
[1,]  0,02500233  0,02648121 -0,06232153
[2,]  0,02648121  0,21739805 -0,28297823
[3,] -0,06232153 -0,28297823  0,40396359
In [3]:
#==============================================================================
# Aproveitamos o exemplo para verificar:
#==============================================================================
print("==== Proporção da variabilidade ====")
print(delta/3)
print("==== Proporção da variabilidade acumulada ====")
print(cumsum(delta)/3)

[1] "==== Proporção da variabilidade ===="
[1] 0,78454534 0,20533887 0,01011579
[1] "==== Proporção da variabilidade acumulada ===="
[1] 0,7845453 0,9898842 1,0000000
In [4]:
#==============================================================================
# Para conferir, note que delta_k = ell.1_k'*ell.1_k
# Logo, podemos ver a proporcao da variabilidade devida ao k-esimo fator 
# atraves da soma de quadrados dos valores da k-esima coluna de L
#==============================================================================
print(delta)
print(cbind((ell.1)%*%ell.1,(ell.2)%*%ell.2,(ell.3)%*%ell.3))

[1] 2,35363603 0,61601660 0,03034736
         [,1]      [,2]       [,3]
[1,] 2,353636 0,6160166 0,03034736