aula de hoje:

1 estimativa de densidades: histogramas

2 estimativa de densidades: KDE

3 análise de agrupamento: K-means

4 PCA: análise de componentes principais

5 LLE: locally linear embbeding

6 o pacote dimRed e algoritmos de redução de dimensionalidade

Para saber o tempo de processamento deste Rmd:

# tempo de processamento do Rmd - inicialização:
t00 = Sys.time()

1. estimativa de densidades: histogramas

temos um conjunto de N dados que queremos representar com um histograma: qual o tamanho ``ideal’’ do bin?

regra de Scott: \(\sigma_s\) - dispersão da amostra \[ \Delta = \frac{3.5 \sigma_s}{N^{1/3}}\]

regra de Freedman-Diaconis: \[ \Delta = \frac{2(q_{75} - q_{25})}{N^{1/3}} = \frac{2.7 \sigma_G}{N^{1/3}} \]

regra de Silverman: \[ \Delta = \frac{1.06 \sigma_s}{N^{1/5}} \]

regra de Knuth: (https://arxiv.org/pdf/physics/0605197.pdf) - método bayesiano - posterior do número de bins M com contagens \(n_k\) em cada bin: \[ P(M|D) = N \log M + \log \Bigg[ \Gamma \Bigg(\frac{M}{2} \Bigg) \Bigg] - M \log \Bigg[ \Gamma \Bigg(\frac{1}{2} \Bigg) \Bigg] - \] \[ -\log \Bigg[ \Gamma \Bigg(\Gamma(N+\frac{M}{2} \Bigg) \Bigg] + \sum_{k=1}^M \log \Bigg[ \Gamma \Bigg(n_k+\frac{1}{2} \Bigg) \Bigg] \] \(M\) é obtido maximizando-se o posterior

Vamos considerar agora uma amostra de elipticidades \((1-b/a)\) de galáxias de baixo brilho superficial na amostra do DR1 do S-PLUS (Mendes de Oliveira et al. 2019):

# distribuições normalizadas:
tabela = read.table(file="ell.dat", header=TRUE)

# dimensões do arquivo:
dim(tabela)
## [1] 244   1
head(tabela)
##   ellipticity
## 1  0.09873418
## 2  0.59826947
## 3  0.52577320
## 4  0.27760252
## 5  0.15665796
## 6  0.40070922
ell = tabela$ellipticity
par(mfrow = c(1,2))
hist(ell,xlab = '1-B/A', col='red',main='',breaks=40)

# distribuições normalizadas: histograma com freq = FALSE
hist(ell,xlab = '1-B/A', col='red',main='',breaks=40,freq=FALSE)

tamanho dos bins:

# regra de Scott:
sigma = sd(ell)
bin_scott = 3.5*sigma/length(ell)^(1/3)
bin_scott
## [1] 0.09861562
# número de bins
ldata = max(ell)-min(ell)
n_scott=ldata/bin_scott+1
# para arredondar:
round(n_scott)
## [1] 8
# regra de Freedman-Diaconis
sig_G = 0.7413*(quantile(ell,0.75,names = FALSE) - quantile(ell,0.25,names = FALSE))
bin_fd = 3.5*sig_G/length(ell)^(1/3)
bin_fd
## [1] 0.115755
n_fd = ldata/bin_fd+1
round(n_fd)
## [1] 7
# regra de Silverman:
bin_silverman = 1.06*sigma/length(ell)^(1/5)
bin_silverman
## [1] 0.06215873
n_silverman = ldata/bin_silverman+1
round(n_silverman)
## [1] 13
# regra de Knuth: teste com até 100 bins

# regra de Knuth
# https://arxiv.org/pdf/physics/0605197.pdf
# optBINS finds the optimal number of bins for a one-dimensionaldata set using the posterior probability for the number of bins
# This algorithm uses a brute-force search trying every possible bin number in the given range.  This can of course be improved.
# Generalization to multidimensional data sets is straightforward.
## Usage:
#           optBINS = function(data,maxM)
# Where:
#           data is a (1,N) vector of data points
#           maxM is the maximum number of bins to consider
## Ref: K.H. Knuth. 2012. Optimal data-based binning for histograms
# and histogram-based probability density models, Entropy.
optBINS = function(data,maxM){

N = length(data)

# Simply loop through the different numbers of bins
# and compute the posterior probability for each.

logp = rep(0,N)

for(M in 2:maxM){

# Bin the data (equal width bins here)
n = table(cut(data, breaks = M))

soma = 0
for(k in 1:M){soma = soma + lgamma(n[k]+0.5)}

logp[M] = N*log(M) + lgamma(M/2) - lgamma(N+M/2) - M*lgamma(1/2)+ soma
#sum(lgamma(n+0.5))
#print(c(M,logp[M]))
}
return(logp)
}

a = optBINS(ell,100)
Mknuth = which(a == max(a))
Mknuth
## [1] 8
# número de bins:
c(round(n_scott),round(n_fd),round(n_silverman),round(Mknuth))
## [1]  8  7 13  8
# visualização
par(mfrow = c(2,2))
plot(cut(ell, breaks = n_scott),col='red')
legend('topleft', 'Scott',cex=1.5,bty='n')
plot(cut(ell, breaks = n_fd),col='green')
legend('topleft', 'Freedman-Diaconis',cex=1.5,bty='n')
plot(cut(ell, breaks = n_silverman),col='darkgoldenrod')
legend('topleft', 'Silverman',cex=1.5,bty='n')
plot(cut(ell, breaks = Mknuth),col='blue')
legend('topleft', 'Knuth',cex=1.5,bty='n')

note que a escolha do tamanho dos bins não é óbvia pois o resultado depende da forma da distribuição!

2. estimativa de densidades: KDE (kernel density estimation)

Kernel: uma função que é aplicada aos dados para modelar/suavizar sua distribuição

método que aplica-se em espaços de dados de qualquer dimensão \(D\)

densidade (pdf) num ponto \(x\): \[ {\hat f_N (x)} = \frac{1}{N h^D} \sum_{i=1}^N K \Bigg( \frac{d(x,x_i)}{h} \Bigg)\]

  • \(K(u)\): função kernel
  • \(d(x_1,x_2)\): ``distância’’ entre \(x_1\) e \(x_2\)
  • \(h\): largura de banda - define o tamanho do kernel

\(h\): não pode ser muito pequeno (a variância seria muito grande), nem muito grande (o viés seria muito grande)

parâmetros livres do KDE: tipo de kernel; largura do kernel

tipos de kernel:

-kernel gaussiano: \[ K(u) = \frac{1}{(2 \pi)^{D/2}}e^{-u^2/2}, ~~~~~~~~ u = d(x,x_i)/h\]

-kernel de Epanechnikov: ótimo em termos da mínima variança \[K(u) = \begin{cases} \frac{3}{4}(1-u^2) & \quad \text{se } |u| \le 1 \\ 0 & \quad \text{se } |u| > 1 \end{cases} \]

-kernel cartola (top-hat): \[K(u) = \begin{cases} \frac{1}{V_D(1)} & \quad \text{se } u \le 1 \\ 0 & \quad \text{se } u > 1 \end{cases} \]

\(V_D(r)\)- volume de uma hiper-esfera \(D\)-dimensional de raio \(r\): \[ V_D(r) = \frac{2 \pi^{D/2} r^D}{D \Gamma(D/2)} \]

Vamos ver os tipos de kernel disponíveis no R:

require(graphics)

# tipos de kernels:
kernels <- eval(formals(density.default)$kernel)
kernels 
## [1] "gaussian"     "epanechnikov" "rectangular"  "triangular"   "biweight"    
## [6] "cosine"       "optcosine"

forma dos kernels:

par(mfrow = c(1,1))
plot (density(0, bw = 1), xlab = "", main = "kernels com bw = 1")
for(i in 2:length(kernels)){
    lines(density(0, bw = 1, kernel =  kernels[i]), col = i)
    legend(1.5,.4, legend = kernels, col = seq(kernels), lty = 1, cex = .8, y.intersp = 1)}

# para saber mais sobre a função density:
?density

Vamos agora examinar a distribuição de elipticidades das galáxias de baixo brilho superficial.

as figuras abaixo têm valores diferentes de h:

# D=1: aplicação a ell.dat
# efeito da largura de banda
par(mfrow = c(2,2))
h=c(1/80,1/40,1/20,bw.SJ(ell))
## SJ is a sensible automatic choice
for(j in 1:4){
plot(density(ell, bw = h[j]), main = "mesmo h, 7 kernels diferentes")
for(i in 2:length(kernels)){lines(density(ell, bw = h[j], kernel = kernels[i]), col = i)}}

vamos agora determinar a largura de banda por validação cruzada, deixando um de fora.

verossimilhança deixando um de fora: \[ VC(h) = \frac{1}{N} \sum_{i=1}^N \log {\hat f}_{h,-i}(x_i) \] a densidade é estimada deixando-se de lado o i-ésimo dado e \(h\) é obtido minimizando-se \(VC(h)\)

# https://www.r-bloggers.com/cross-validation-for-kernel-density-estimation/
# o default é o kernel gaussiano
X = ell
J=function(h){
  fhat=Vectorize(function(x) density(X,from=x,to=x,n=1,bw=h)$y)
  fhati=Vectorize(function(i) density(X[-i],from=X[i],to=X[i],n=1,bw=h)$y)
  F=fhati(1:length(X))
  return(integrate(function(x) fhat(x)^2,-Inf,Inf)$value-2*mean(F))
}

vx=seq(.01,1,by=.01)
vy=Vectorize(J)(vx)
df=data.frame(vx,vy)

par(mfrow = c(2,2))
library(ggplot2)
qplot(vx,vy,geom="line",data=df)

# optimal value:
a = optimize(J,interval=c(.01,1))
a
## $minimum
## [1] 0.07017881
## 
## $objective
## [1] -1.492248
hopt = a$minimum
hopt
## [1] 0.07017881
hopt=round(hopt,3)

# histograma com o valor ótimo de h:
par(mfrow = c(1,2))
hist(ell,xlab = '1-B/A', col='red',main='',breaks=floor(2/hopt))
y = density(ell, bw = hopt)
plot(y,xlab = '1-B/A', col='red', type = "n", main = sprintf(paste('h = ',hopt)))
polygon(y, col = "wheat")

Vamos agora considerar o KDE em duas dimensões:

Vamos examinar a distribuição das cores (g-r) x (r-i) de uma amostra de galáxias extraída do S-PLUS.

# kde2d: two-dimensional kernel density estimation with an axis-aligned bivariate normal kernel, evaluated on a square grid

library(MASS)
library(KernSmooth)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
tabela = read.table(file="rmixgmr.dat", header=TRUE)
dim(tabela)
## [1] 1025    2
head(tabela)
##    rmi  gmr
## 1 0.28 0.60
## 2 0.18 0.58
## 3 0.05 0.46
## 4 0.19 0.43
## 5 0.19 0.39
## 6 0.38 0.76
rmi = tabela$rmi
gmr = tabela$gmr

# vamos olhar a cara das distribuições:
par(mfrow = c(1,2))
hist(rmi,xlab='r-i',main='',col = 'red')
hist(gmr,xlab='g-r',main='',col = 'red')

# resultados para vários valores de h
par(mfrow = c(2,2))
f <- kde2d(rmi,gmr, n = 500, h = 0.02)
image(f,xlab='r-i',ylab='g-r',main='h = 0.02')
f <- kde2d(rmi,gmr, n = 500, h = 0.025)
image(f,xlab='r-i',ylab='g-r',main='h = 0.025')
f <- kde2d(rmi,gmr, n = 500, h = 0.05)
image(f,xlab='r-i',ylab='g-r',main='h = 0.05')
f <- kde2d(rmi,gmr, n = 500, h = 0.075)
image(f,xlab='r-i',ylab='g-r',main='h = 0.075')

plot(rmi,gmr, pch=20,xlab='r-i',ylab='g-r')
f <- kde2d(rmi,gmr, n = 500, h = 0.1)
image(f,xlab='r-i',ylab='g-r',main='h = 0.1')
f <- kde2d(rmi,gmr, n = 500, h = 0.2)
image(f,xlab='r-i',ylab='g-r',main='h = 0.2')
f <- kde2d(rmi, gmr, n = 500, h = c(width.SJ(rmi), width.SJ(gmr)) )
image(f,xlab='r-i',ylab='g-r',main='h = SJ')

# a largura de banda S-J (Sheather e Jones 1991) é frequentemente considerada a "melhor escolha"
c(width.SJ(rmi), width.SJ(gmr)) 
## [1] 0.1074551 0.1615408

3. análise de agrupamento: K-means

tipos de agrupamento:

-``duro’’: cada objeto pertence ou não a um dado grupo

-``suave’’: cada objeto pertence a um dado grupo com uma certa probabilidade

exemplos de métodos:

-algoritmos combinatórios: K-means, hierarchical clustering

-modelos probabilísticos: misturas gaussianas,…

tipos de distâncias entre 2 objetos (\(\mathbf{x}_i\) e \(\mathbf{x}_j\)):

-euclidiana: \(d_{ij} = \sum_{k=1}^M |\mathbf{x}_i-\mathbf{x}_j|\)

-manhattan: \(d_{ij} = \Big[\sum_{k=1}^M |\mathbf{x}_i-\mathbf{x}_j|^2\Big]^{1/2}\)

-mahalanobis: $d_{ij} = ^{1/2} $

algoritmo K-means:

-objetivo: atribuir N objetos a K grupos, K\(< <\)N \(\longrightarrow\) K é dado!

o algoritmo procura, iterativamente, minimizar a soma \[ J = \sum_{i=1}^N \sum_{k=1}^K (\mathbf{x}_{i}-\mathbf{\mu}_{k})^2 = \sum_{i=1}^N \sum_{j=1}^M \sum_{k=1}^K (x_{ij}-\mu_{jk})^2\] onde \(\mathbf{\mu}_k\) é o centróide do grupo K

algoritmo:

-1. defina o número de grupos, K

-2. inicialize, escolhendo K objetos ao acaso como centro

-3. atribua cada objeto ao grupo mais próximo e recalcule o centróide do grupo

repita 3 até a convergência (os grupos ‘congelam’)

K deve ser dado a priori

o algoritmo é sensível a inicialização: diferentes ‘rodadas’ podem levar a grupos diferentes

Vamos considerar inicialmente dados simulados em duas dimensões:

# vamos gerar coordenadas (x,y) para dois conjuntos de dados
set.seed(123)
# grupo 1: mu = (-1,-1)    sd = 0.5 
mu1 = c(-1,-1)
sig = 0.5
grupo1 = cbind(rnorm(50,mu1[1],sig), rnorm(50,mu1[2],sig))
# grupo 2: mu = (-1,-1)    sd = 0.5 
mu2 = c(1,1)
sig = 0.5
grupo2 = cbind(rnorm(50,mu2[1],sig), rnorm(50,mu2[2],sig))

dados = rbind(grupo1, grupo2)

# a função kmeans do R usa um algoritmo modificado que reinicializa 
# aleatoriamente o algoritmo várias vezes
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos
kmres$centers[1,]
## [1] -0.9827982 -0.9267959
kmres$centers[2,]
## [1] 0.8730498 1.0194034
#visualização
par(mfrow = c(1,1))
plot(dados, pch = 16, xlab = "x", ylab = "y", col = kmres$cluster,xlim=c(-5,5),ylim=c(-5,5))
points(kmres$centers[,1], kmres$centers[,2],pch=4,lwd=4, col="blue", cex = 1.5)
points(mu1[1],mu1[2],pch=4,lwd=4, col="green", cex = 1.5)
points(mu2[1],mu2[2],pch=4,lwd=4, col="green", cex = 1.5)
legend('topleft', col=c("green", 'blue'), lwd=2, legend=c("input", "k-means"))

# número de objetos em cada grupo:
kmres$size
## [1] 50 50
# classificação de cada objeto:
kmres$cluster
##   [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##  [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##  [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
# vamos ver como fica com sd = 1.5 
sig = 1.5
grupo1 = cbind(rnorm(50,mu1[1],sig), rnorm(50,mu1[2],sig))
grupo2 = cbind(rnorm(50,mu2[1],sig), rnorm(50,mu2[2],sig))
# dados
dados = rbind(grupo1, grupo2)

# função kmeans
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos
kmres$centers[1,]
## [1] -1.5149972 -0.5467696
kmres$centers[2,]
## [1] 1.4648223 0.4836374
# plots
plot(dados, pch = 16, xlab = "x", ylab = "y", col = kmres$cluster,xlim=c(-5,5),ylim=c(-5,5))
points(kmres$centers[,1], kmres$centers[,2],pch=4,lwd=4, col="blue", cex = 1.5)
points(mu1[1],mu1[2],pch=4,lwd=4, col="green", cex = 1.5)
points(mu2[1],mu2[2],pch=4,lwd=4, col="green", cex = 1.5)
legend('topleft', col=c("green", 'blue'), lwd=2, legend=c("input", "k-means"))

# número de objetos em cada grupo:
kmres$size
## [1] 47 53
# classificação de cada objeto:
kmres$cluster
##   [1] 1 1 1 2 1 2 1 1 1 1 1 1 2 1 1 1 1 2 1 1 2 1 1 1 2 1 1 1 1 2 1 1 1 1 1 1 1
##  [38] 1 2 2 2 1 1 1 1 2 2 1 1 1 2 2 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 1 2 2 1
##  [75] 2 2 1 2 2 2 2 2 2 1 2 2 2 2 2 2 1 2 2 2 2 2 1 2 2 2

Vamos usar a rotina kernel2d, que ajusta um kernel em 2 dimensões:

# visualização dos dados com kernel2d
par(mfrow = c(1,1))
# Two-dimensional kernel-density estimator
library(MASS)
dados_sm = kde2d(rmi, gmr, h=0.075,n=500)
image(dados_sm, col=grey(13:0/15), xlab='r-i',ylab='g-r',main='galáxias do S-PLUS')

# padronização das variáveis
dados=cbind(rmi,gmr)
# o default do scale é: scale(x, center = TRUE, scale = TRUE)
dados_scale = scale(dados)

# K-means variando o número de grupos entre 2 e 15
# o atributo withinss do K-means é o vetor com a soma dos quadrados das distâncias dentro de cada grupo
wss = (nrow(dados_scale)-1)*sum(apply(dados_scale,2,var)) 
for (i in 2:15){ wss[i] = sum(kmeans(dados_scale,centers=i)$withinss)}

# visualização
par(mfrow=c(1,1))
plot(1:15, wss, type="b", xlab="número de grupos", ylab="soma dos quadrados das distâncias dentro de cada grupo") 

# ajuste de 2 grupos
kmres = kmeans(dados_scale, 2)

# número de objetos em cada grupo:
kmres$size 
## [1] 560 465
# médias dos grupos
a = aggregate(dados,by=list(kmres$cluster),FUN=mean)
xm=c(a[1,2],a[1,3])
ym=c(a[2,2],a[2,3])

# adiciono a cada objeto seu grupo 
dados = data.frame(dados, kmres$cluster)
tipos =dados[,3]

#visualização dos grupos
par(mfrow = c(1,2))
# Two-dimensional kernel-density estimator
dados_sm = kde2d(rmi, gmr, h=0.075,n=500)
image(dados_sm, col=grey(13:0/15), xlab='r-i',ylab='g-r',main='')
points(xm[1],xm[2],pch=4,lwd=4, col="green", cex = 1.5)
points(ym[1],ym[2],pch=4,lwd=4, col="green", cex = 1.5)

image(dados_sm, col=grey(13:0/15), xlab='r-i',ylab='g-r',main='')
points(dados[tipos == 1,1], dados[tipos == 1,2],col='red')
points(dados[tipos == 2,1], dados[tipos == 2,2],col='blue')
points(xm[1],xm[2],pch=4,lwd=4, col="green", cex = 1.5)
points(ym[1],ym[2],pch=4,lwd=4, col="green", cex = 1.5)

4. PCA - Análise de Componentes Principais

para que serve?

-compressão de dados, via redução da dimensionalidade

-visualização de dados multidimensionais

como funciona?

  1. no espaço dos dados identifique a direção que apresenta a maior variância dos dados: esta é a direção da primeira componente principal, \(PC_1\)
  1. no sub-espaço perpendicular a \(PC_1\) identifique a direção que apresenta a maior variância dos dados: esta é a direção da segunda componente principal, \(PC_2\)
  1. no sub-espaço perpendicular a \(PC_1\) e \(PC_2\) identifique…
  1. repita o procedimento até \(PC_D\)

PCA cria um novo sistema de coordenadas no espaço de dados

em geral poucas componentes são responsáveis pela maior parte da variância

frequentemente os objetos podem ser descritos por um número \(d << D\) de componentes

como o PCA funciona?

vamos considerar um conjunto de N objetos, cada um descrito por D variáveis:

\(x_{ij}\): variável \(j\) do \(i\)-ésimo objeto

matriz de covariância dos dados: matriz \(D \times D\) definida como \[ C_{jk} = \frac{1}{N} \sum_{i=1}^N (x_{ij}-{\bar x}_j) (x_{ik}-{\bar x}_k) \] onde \[ {\bar x}_j = \frac{1}{N} \sum_{i=1}^N x_{ij} \] note que as variâncias dos dados correspondem à diagonal, \(C_{jj}\)

muitas vezes as medidas das variáveis são incompatíveis (como distâncias e magnitudes):

para ``resolver’’ isso padronizam-se as variáveis, convertendo-as em quantidades de média nula e variância unitária: \[ x_{ij}^\prime = \frac{(x_{ij}-{\bar x}_j)}{\sqrt{C_{jj}}} \] e constrói-se a matriz de covariância com as variáveis padronizadas

cada componente principal é descrita como uma combinação linear das variáveis padronizadas: \[ PC_j = \sum_{i=1}^N \sum_{k=1}^D a_{jki} x_{ik}^\prime = \sum_{i=1}^N \mathbf{a}_{ji} \mathbf{x_i^\prime} \]

os \(\mathbf{a}_j\) são os autovetores de \(C\)

o autovetor de \(C\) com o j-ésimo maior autovalor é a j-ésima componente principal e satisfaz (como as demais componentes):
\[\mathbf{a}_j\mathbf{a}_k^\prime = \delta_{jk}\] (normalização e ortogonalidade; \(\mathbf{a}_k^\prime\) é o vetor transposto)

os autovalores são proporcionais à variância de cada componente

assim, para descrever os dados pode-se selecionar o número de componentes em termos da fração da variância total que se deseja ``explicar’’

problemas com o PCA

  • as PC não têm necessariamente significado físico
  • a seleção do número de componentes de interesse é muitas vezes subjetiva
  • a análise depende do tipo de normalização aplicada aos dados
  • é um método linear

Vamos ilustrar o PCA inicialmente com dados do levantamento fotométrico S-PLUS, que utiliza um sistema de 12 filtros para obter bons redshifts fotométricos e estudar evolução de galáxias. Aqui vamos considerar apenas as cores (11) de galáxias com redshifts fotométricos entre 0.10 e 0.15 com fotometria em todas as bandas. A precisão dos redshifts fotométricos do S-PLUS é \(\sim \Delta_z/(1+z) = 0.02\)

inicialização

set.seed(123)
library(lle)
## Loading required package: scatterplot3d
## Loading required package: snowfall
## Loading required package: snow
library(plot3D)
library(KernSmooth)
library(dimRed)
## Loading required package: DRR
## Loading required package: kernlab
## 
## Attaching package: 'kernlab'
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## Loading required package: CVST
## Loading required package: Matrix
## 
## Attaching package: 'dimRed'
## The following object is masked from 'package:stats':
## 
##     embed
## The following object is masked from 'package:base':
## 
##     as.data.frame
library(fastICA)
library(RSpectra)
library(corrplot)
## corrplot 0.84 loaded
library(reticulate)

t0 = Sys.time()

leitura dos dados do S-PLUS: cores relativas à banda r_petro

cores = read.table('cores010015.dat',header=T, fill=T)

Algumas características destes dados:

dim(cores)
## [1] 1575   11
head(cores)
##   uJAVA_petro.r_petro F378_petro.r_petro F395_petro.r_petro F410_petro.r_petro
## 1                2.20               1.74               1.25               1.34
## 2                1.17               2.22               1.02               0.57
## 3                1.23               0.94               0.92               0.69
## 4                1.96               1.63               1.71               1.37
## 5                1.82               1.61               0.93               0.36
## 6                2.00               2.04               2.01               0.91
##   F430_petro.r_petro g_petro.r_petro F515_petro.r_petro F660_petro.r_petro
## 1               0.65            0.56               0.26              -0.32
## 2               0.54            0.31               0.22              -0.09
## 3               0.31            0.27               0.49               0.01
## 4               0.97            0.67               0.33              -0.16
## 5               0.74            0.32               0.48              -0.09
## 6               0.89            0.71               0.26              -0.19
##   i_petro.r_petro F861_petro.r_petro z_petro.r_petro
## 1           -0.38              -0.65           -0.74
## 2           -0.24              -0.41           -0.47
## 3            0.16              -0.28           -0.07
## 4           -0.35              -0.49           -0.56
## 5           -0.23              -0.35           -0.42
## 6           -0.36              -0.62           -0.50

Vamos agora calcular o PCA, padronizando as variáveis com média zero e variância unitária:

pca.cores = prcomp(cores, retx=TRUE, center=TRUE, scale=TRUE)

# raiz quadrada dos autovalores
sd = pca.cores$sdev

# autovalores
lambda = sd^2

# fraçao da variância explicada
round(cumsum(lambda)/sum(lambda),digits = 2)
##  [1] 0.38 0.55 0.62 0.69 0.75 0.80 0.85 0.90 0.94 0.97 1.00
# autovetores
autovec = pca.cores$rotation

# componentes principais
pc = pca.cores$x

# visualização: scree plot
screeplot(pca.cores, main=" scree plot")

# visualização: duas primeiras componentes
par(mfrow = c(1,1))
plot(pc[,1], pc[,2], xlab="PCA 1", ylab="PCA 2",asp=1,pch=20,col='red')

par(mfrow = c(2,2))
smoothScatter(pc[,1],pc[,2],nrpoints=0,add=FALSE, xlab='pc1', ylab='pc2', main = '',asp = 1)

smoothScatter(pc[,1],pc[,3],nrpoints=0,add=FALSE, xlab='pc1', ylab='pc3', main = '',asp = 1)

smoothScatter(pc[,1],pc[,4],nrpoints=0,add=FALSE, xlab='pc1', ylab='pc4', main = '',asp = 1)

smoothScatter(pc[,2],pc[,3],nrpoints=0,add=FALSE, xlab='pc2', ylab='pc3', main = '',asp = 1)

distância de Mahalanobis: distância entre dois pontos num espaço multivariado que leva em conta a correlação entre variáveis: \[ D_{M}(x)={\sqrt {(x-\mu )^{T}C^{-1}(x-\mu )}} \] \(C\): matriz de covariância

ndim=dim(cores)[2]
for (j in 1:ndim){
media= mean(cores[,j])
sdj=sd(cores[,j])
mah = mahalanobis(cores,center=T,cov=cov(cores))}

par(mfrow = c(1,1))
id = seq(1,dim(cores)[1],1)
lmah = log10(mah)
plot(id,lmah,main = 'log da distância de  mahalanobis', xlab= 'obj id', ylab = 'log(D)')

esta distância muitas vezes é usada para se identificar outliers

vamos agora examinar a relação entre as cores das galáxias e as duas primeiras componentes principais via análise de correlação:

# primeiro criamos uma matriz com as cores e as duas componentes principais do PCA:
aux = cbind(pc[,1],pc[,2],cores)
# análise de correlação
correlacao = cor(aux, method = "spearman")
#jpeg('corr.jpg')
par(mfrow = c(1,1))
corrplot(correlacao, type = "upper",  tl.col = "black", tl.srt = 45)

# coeficientes de correlação:
aux = as.matrix(cbind(pc[,1],pc[,2],cores))
aux = cbind(pc[,1],pc[,2],cores)
correlacao = cor(aux, method = "spearman")
#jpeg('corr.jpg')
par(mfrow = c(1,1))
corrplot(correlacao, type = "upper",  tl.col = "black", tl.srt = 45)

> Note que as duas primeiras linhas na figura acima mostram justamente as correlações/anti-correlações das PCs 1 e 2 com as cores.

5. LLE: locally linear embbeding

embedding: : quando uma estrutura está contida dentro de outra

representa dados de altas dimensões em um espaço de dimensão menor, preservando a vizinhança local de cada ponto

vizinhança: k-vizinhos mais próximos

algoritmo não-supervisionado

o algoritmo (iterativo) tem 2 passos:

  1. para cada ponto determina-se um conjunto de pesos \(W\) que melhor o reconstroi a partir de seus k-vizinhos mais próximos (ajuste de um hiperplano)- esses pesos codificam a geometria local
  1. com os pesos fixos, determina-se um conjunto de dados de menor dimensão que preserva as relações de vizinhança desses pesos

LLE: \(W\) obtido via álgebra linear

vamos considerar inicialmente a curva S em 3 dimensões:

data(lle_scurve_data)
X <- lle_scurve_data
head(X, n=6)
##            [,1]      [,2]       [,3]
## [1,]  0.9553277 4.9512047 -0.1744640
## [2,] -0.6600097 3.2669075 -0.7732213
## [3,] -0.9832013 1.2573407 -0.2963627
## [4,]  0.9543690 1.6825490 -0.1797167
## [5,]  0.9576142 0.1857008 -0.1612359
## [6,]  0.8520630 0.5582119 -0.4707255
# plot 3-D dos dados
par(mfrow = c(1,1))
scatterplot3d(x=X[,1], y=X[,2], z=X[,3], color=X[,2])

a análise depende do número de vizinhos de cada ponto

se o número de vizinhos é pequeno, o mapeamento não vai refletir nenhuma propriedade global, enquanto que se for muito grande vais perder seu caráter não-linear e se comportar como PCA

vamos inicialmente ver o resultado para vários valores do número de vizinhos k:

# K=5
# m: dimensão intrínseca dos dados- para visualização; a dimensão intrínseca é calculada pelo programa
# reg: método de regularização
# v: 
k5 <- lle(X, m=2, k=5, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=1.92125, mode=2
## computing coordinates
plot(k5$Y, main="K=5 data", xlab=expression(y[1]), ylab=expression(y[2]))

# K=15
k15 <- lle(X, m=2, k=15, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=2, mode=2
## computing coordinates
plot(k15$Y, main="K=15 data", xlab=expression(y[1]), ylab=expression(y[2]))

# K=40
k40 <- lle(X, m=2, k=40, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=2.00625, mode=2
## computing coordinates
plot(k40$Y, main="K=40 data", xlab=expression(y[1]), ylab=expression(y[2]))

o pacote lle tem uma rotina para escolher o “melhor” valor de k:

# deve-se escolher k que dá o menor valor de 1-rho^2
# notem que estou usando as 8 CPUs de meu laptop!
par(mfrow = c(1,1))
calc_k(X, m=2, kmin=1, kmax=40, plotres=TRUE, parallel=TRUE, cpus=8, iLLE=FALSE)
## Warning in searchCommandline(parallel, cpus = cpus, type = type, socketHosts =
## socketHosts, : Unknown option on commandline: rmarkdown::render('/home/laerte/
## cursos/dados/2020A/aula11/R/aula11.Rmd',~+~~+~encoding~+~
## R Version:  R version 3.6.3 (2020-02-29)
## snowfall 1.84-6.1 initialized (using snow 0.4-3): parallel execution on 8 CPUs.
## Library lle loaded.
## Library lle loaded in cluster.
## 
## Stopping cluster

##     k       rho
## 1   1 0.9987824
## 2   2 0.9876860
## 3   3 0.9925197
## 4   4 0.9127051
## 5   5 0.6969111
## 6   6 0.6765670
## 7   7 0.4045646
## 8   8 0.4181170
## 9   9 0.3045077
## 10 10 0.2730377
## 11 11 0.1722015
## 12 12 0.1522597
## 13 13 0.1589914
## 14 14 0.1489428
## 15 15 0.1605777
## 16 16 0.1512468
## 17 17 0.1468803
## 18 18 0.2965946
## 19 19 0.1507325
## 20 20 0.1694365
## 21 21 0.2238161
## 22 22 0.1745816
## 23 23 0.2007847
## 24 24 0.2011773
## 25 25 0.2389543
## 26 26 0.2231868
## 27 27 0.3097121
## 28 28 0.2208020
## 29 29 0.4766895
## 30 30 0.2791912
## 31 31 0.1252955
## 32 32 0.1472460
## 33 33 0.1350042
## 34 34 0.1492596
## 35 35 0.1546612
## 36 36 0.1470390
## 37 37 0.1288793
## 38 38 0.1175156
## 39 39 0.1228238
## 40 40 0.1176667

escolhendo k=31:

k31 <- lle(X, m=2, k=31, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=2, mode=2
## computing coordinates
plot(k31$Y, main="K=31 data", xlab=expression(y[1]), ylab=expression(y[2]))

uma das saídas do lle é o “número intrínseco de dimensões”: isso é definido pelo número de componentes que “explicam” uma fração \(v\) da variância; em geral se usa \(v\) entre 0.9 e 0.99

Vamos agora analisar nossa amostra com as cores das galáxias do S-PLUS. Vamos começar vendo os resultados para alguns valores de k:

k5 <- lle(cores, m=2, k=5, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=3.364444, mode=3
## computing coordinates
plot(k5$Y, main="K=5 data", xlab=expression(y[1]), ylab=expression(y[2]))

# K=15
k15 <- lle(cores, m=2, k=15, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=5.069841, mode=5
## computing coordinates
plot(k15$Y, main="K=15 data", xlab=expression(y[1]), ylab=expression(y[2]))

# K=40
k40 <- lle(cores, m=2, k=40, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=5.653333, mode=6
## computing coordinates
plot(k40$Y, main="K=40 data", xlab=expression(y[1]), ylab=expression(y[2]))

qual seria o melhor valor de k?

calc_k(cores, m=2, kmin=1, kmax=40, plotres=TRUE, parallel=TRUE, cpus=8, iLLE=FALSE)
## Warning in searchCommandline(parallel, cpus = cpus, type = type, socketHosts =
## socketHosts, : Unknown option on commandline: rmarkdown::render('/home/laerte/
## cursos/dados/2020A/aula11/R/aula11.Rmd',~+~~+~encoding~+~
## snowfall 1.84-6.1 initialized (using snow 0.4-3): parallel execution on 8 CPUs.
## Library lle loaded.
## Library lle loaded in cluster.
## 
## Stopping cluster

##     k       rho
## 1   1 0.9980250
## 2   2 0.9877269
## 3   3 0.8519231
## 4   4 0.8395129
## 5   5 0.7685058
## 6   6 0.8027322
## 7   7 0.7554615
## 8   8 0.8276032
## 9   9 0.8679262
## 10 10 0.8752581
## 11 11 0.7889229
## 12 12 0.8468184
## 13 13 0.7057745
## 14 14 0.7284917
## 15 15 0.6521474
## 16 16 0.5215134
## 17 17 0.4460243
## 18 18 0.4458253
## 19 19 0.5474950
## 20 20 0.4220398
## 21 21 0.3506301
## 22 22 0.4803771
## 23 23 0.3790212
## 24 24 0.3218850
## 25 25 0.3855203
## 26 26 0.3637025
## 27 27 0.4670727
## 28 28 0.3751212
## 29 29 0.4447791
## 30 30 0.3710600
## 31 31 0.3207646
## 32 32 0.3997840
## 33 33 0.3437415
## 34 34 0.3808011
## 35 35 0.3817986
## 36 36 0.4371570
## 37 37 0.4006051
## 38 38 0.4109186
## 39 39 0.3420062
## 40 40 0.3489851

para k=24:

k24 <- lle(cores, m=2, k=24, reg=2, ss=FALSE, id=TRUE, v=0.9 )
## finding neighbours
## calculating weights
## intrinsic dim: mean=5.466032, mode=6
## computing coordinates
plot(k24$Y, main="K=24 data", xlab=expression(y[1]), ylab=expression(y[2]))

plot(k24$id, main="intrinsic dimension", type="l", xlab=expression(x[i]), ylab="id", lwd=2 )

vamos agora examinar a correlação entre as cores das galáxias e suas coordenadas no espaço reduzido:

# criamos uma matriz com as cores e as duas dimensões do espaço projetado:
aux = cbind(k24$Y,cores)
# análise de correlação
correlacao = cor(aux, method = "spearman")
par(mfrow = c(1,1))
corrplot(correlacao, type = "upper",  tl.col = "black", tl.srt = 45)

correlacao
##                               1           2 uJAVA_petro.r_petro
## 1                    1.00000000  0.09382637           0.6769820
## 2                    0.09382637  1.00000000           0.3616589
## uJAVA_petro.r_petro  0.67698204  0.36165887           1.0000000
## F378_petro.r_petro   0.25315451  0.75734277           0.3896224
## F395_petro.r_petro   0.82548166  0.28296875           0.3649365
## F410_petro.r_petro   0.18784748  0.70597318           0.3757954
## F430_petro.r_petro   0.34291068  0.77402692           0.4483620
## g_petro.r_petro      0.42703354  0.63967239           0.5444221
## F515_petro.r_petro   0.34975798  0.40953160           0.3606841
## F660_petro.r_petro  -0.08291645 -0.12263624          -0.1067624
## i_petro.r_petro     -0.26122251 -0.41973686          -0.3223626
## F861_petro.r_petro  -0.28855249 -0.43682312          -0.3567251
## z_petro.r_petro     -0.27476398 -0.51628293          -0.3719825
##                     F378_petro.r_petro F395_petro.r_petro F410_petro.r_petro
## 1                            0.2531545         0.82548166         0.18784748
## 2                            0.7573428         0.28296875         0.70597318
## uJAVA_petro.r_petro          0.3896224         0.36493647         0.37579541
## F378_petro.r_petro           1.0000000         0.36150573         0.41782128
## F395_petro.r_petro           0.3615057         1.00000000         0.35728748
## F410_petro.r_petro           0.4178213         0.35728748         1.00000000
## F430_petro.r_petro           0.4616368         0.41635390         0.47561127
## g_petro.r_petro              0.4939474         0.43303139         0.53214639
## F515_petro.r_petro           0.3580908         0.33037524         0.32081262
## F660_petro.r_petro          -0.0347785        -0.08794922        -0.09885729
## i_petro.r_petro             -0.2805226        -0.28910284        -0.32850582
## F861_petro.r_petro          -0.2917768        -0.29126457        -0.31499528
## z_petro.r_petro             -0.3453992        -0.31870845        -0.39542018
##                     F430_petro.r_petro g_petro.r_petro F515_petro.r_petro
## 1                           0.34291068      0.42703354          0.3497580
## 2                           0.77402692      0.63967239          0.4095316
## uJAVA_petro.r_petro         0.44836203      0.54442211          0.3606841
## F378_petro.r_petro          0.46163682      0.49394739          0.3580908
## F395_petro.r_petro          0.41635390      0.43303139          0.3303752
## F410_petro.r_petro          0.47561127      0.53214639          0.3208126
## F430_petro.r_petro          1.00000000      0.63539454          0.4242650
## g_petro.r_petro             0.63539454      1.00000000          0.4799864
## F515_petro.r_petro          0.42426497      0.47998641          1.0000000
## F660_petro.r_petro         -0.09918295     -0.08867052          0.1053200
## i_petro.r_petro            -0.37118170     -0.37204490         -0.1325445
## F861_petro.r_petro         -0.39030704     -0.47144753         -0.1848162
## z_petro.r_petro            -0.45835752     -0.52193909         -0.2260435
##                     F660_petro.r_petro i_petro.r_petro F861_petro.r_petro
## 1                          -0.08291645      -0.2612225         -0.2885525
## 2                          -0.12263624      -0.4197369         -0.4368231
## uJAVA_petro.r_petro        -0.10676240      -0.3223626         -0.3567251
## F378_petro.r_petro         -0.03477850      -0.2805226         -0.2917768
## F395_petro.r_petro         -0.08794922      -0.2891028         -0.2912646
## F410_petro.r_petro         -0.09885729      -0.3285058         -0.3149953
## F430_petro.r_petro         -0.09918295      -0.3711817         -0.3903070
## g_petro.r_petro            -0.08867052      -0.3720449         -0.4714475
## F515_petro.r_petro          0.10532002      -0.1325445         -0.1848162
## F660_petro.r_petro          1.00000000       0.5178519          0.4548991
## i_petro.r_petro             0.51785191       1.0000000          0.5936204
## F861_petro.r_petro          0.45489911       0.5936204          1.0000000
## z_petro.r_petro             0.43140362       0.6822189          0.6013348
##                     z_petro.r_petro
## 1                        -0.2747640
## 2                        -0.5162829
## uJAVA_petro.r_petro      -0.3719825
## F378_petro.r_petro       -0.3453992
## F395_petro.r_petro       -0.3187085
## F410_petro.r_petro       -0.3954202
## F430_petro.r_petro       -0.4583575
## g_petro.r_petro          -0.5219391
## F515_petro.r_petro       -0.2260435
## F660_petro.r_petro        0.4314036
## i_petro.r_petro           0.6822189
## F861_petro.r_petro        0.6013348
## z_petro.r_petro           1.0000000

Note que as duas primeiras linhas na figura acima mostram justamente as correlações/anti-correlações das 2 projeções com as cores.

Compare as cores “mais explicativas” nesse caso com as obtidas por PCA

6. o pacote dimRed

O R tem um pacote, o dimRed, que reune os principais algoritmos de redução de dimensionalidade, com uma interface comum para análise e plot.

para saber os métodos disponíveis:

dimRedMethodList()
##  [1] "AutoEncoder"         "DiffusionMaps"       "DRR"                
##  [4] "FastICA"             "KamadaKawai"         "DrL"                
##  [7] "FruchtermanReingold" "HLLE"                "Isomap"             
## [10] "kPCA"                "PCA_L1"              "LaplacianEigenmaps" 
## [13] "LLE"                 "MDS"                 "nMDS"               
## [16] "NNMF"                "PCA"                 "tSNE"               
## [19] "UMAP"

o pacote disponibiliza também vários conjuntos de dados para testes:

## lista dos conjuntos de dados disponíveis:
dataSetList()
##  [1] "Swiss Roll"           "Broken Swiss Roll"    "Helix"               
##  [4] "Twin Peaks"           "Sphere"               "Ball"                
##  [7] "FishBowl"             "3D S Curve"           "variable Noise Helix"
## [10] "Iris"                 "Cube"
## para carregar um conjunto de dados:
curva_s <- loadDataSet("3D S Curve")
plot(curva_s, type = "3vars")

neste pacote a função embed é usada para implementar os métodos de redução de dimensionalidade.

vamos fazer um teste com LLE aplicado à curva S usando 45 vizinhos:

emb_LLE <- embed(curva_s, "LLE", knn = 45)
## finding neighbours
## calculating weights
## computing coordinates
print(emb_LLE)
## Method:
## lle 
## Parameters:
## List of 2
##  $ knn : num 45
##  $ ndim: num 2
plot(emb_LLE, type = "2vars")

como avaliar a qualidade da redução de dimensionalidade?

do manual do dimRed:

Rank based criteria

Q_local, Q_global, and mean_R_nx are quality criteria based on the Co-ranking matrix. Q_local and Q_global determine the local/global quality of the embedding, while mean_R_nx determines the quality of the overall embedding. They are parameter free and return a single number. The object must include the original data. The number returns is in the range [0, 1], higher values mean a better local/global embedding.

quality(emb_LLE,'Q_local')
## [1] 0.5001746
quality(emb_LLE,'Q_global')
##        86 
## 0.2520802
quality(emb_LLE,'mean_R_NX')
## [1] 0.492865

vamos agora testar vários métodos:

# PCA
emb_PCA <- embed(curva_s, "PCA")
print(emb_PCA)
## Method:
## PCA 
## Parameters:
## List of 2
##  $ center: logi TRUE
##  $ scale.: logi FALSE
plot(emb_PCA, type = "2vars")

# Isomap:
emb_isomap <- embed(curva_s, "Isomap", knn = 10)
## 2020-06-02 17:26:08: Isomap START
## 2020-06-02 17:26:08: constructing knn graph
## 2020-06-02 17:26:08: calculating geodesic distances
## 2020-06-02 17:26:09: Classical Scaling
print(emb_isomap)
## Method:
## Isomap 
## Parameters:
## List of 4
##  $ knn     : num 10
##  $ ndim    : num 2
##  $ get_geod: logi FALSE
##  $ eps     : num 0
plot(emb_isomap, type = "2vars")

# t-SNE
emb_tsne <- embed(curva_s, "tSNE", perplexity = 80)
print(emb_tsne)
## Method:
## tsne 
## Parameters:
## List of 4
##  $ d         :function (x, method = "euclidean", diag = FALSE, upper = FALSE, p = 2)  
##  $ perplexity: num 80
##  $ theta     : num 0.5
##  $ ndim      : num 2
plot(emb_tsne, type = "2vars")

vamos avaliar a qualidade da redução de dimensionalidade:

quality(emb_PCA,'Q_local')
## [1] 0.513354
quality(emb_LLE,'Q_local')
## [1] 0.5001746
quality(emb_isomap,'Q_local')
## [1] 0.8404466
quality(emb_tsne,'Q_local')
## [1] 0.8285775
# ###
quality(emb_PCA,'Q_global')
##      433 
## 0.328098
quality(emb_LLE,'Q_global')
##        86 
## 0.2520802
quality(emb_isomap,'Q_global')
##        60 
## 0.4020893
quality(emb_tsne,'Q_global')
##        33 
## 0.3827952
# ###
quality(emb_PCA,'mean_R_NX')
## [1] 0.7500491
quality(emb_LLE,'mean_R_NX')
## [1] 0.492865
quality(emb_isomap,'mean_R_NX')
## [1] 0.6753742
quality(emb_tsne,'mean_R_NX')
## [1] 0.6349615
# tempo de processamento do Rmd:
Sys.time() -t00
## Time difference of 5.02563 mins