aula 9

aula de hoje: métodos de aprendizagem de máquina para análise de dados não paramétrica

1 estimativa de densidades:

1.1 KDE kernel density estimation

1.2 GMM gaussian mixture models

2 análise de agrupamento:

2.1 K-means

2.2 GMM

3 PCA: análise de componentes principais

4 redução de dimensionalidade

4.1 PCA

4.2 UMAP: Uniform manifold approximation and projection

5. Exercícios

---

---

---

1. estimativa de densidades

Vamos considerar como modelar uma distribuição em uma e em várias dimensões

Quando normalizadas, as distribuições de densidades podem ser consideradas funções de distribuição de probabilidades

O histograma é o caso mais simples de distribuição de densidade de uma variável aleatória: ele mostra a frequência da variável em intervalos discretos de valores, geralmente constantes, que denominaremos bins

Como vimos, temos nesse caso um problema: qual é o tamanho “ideal” do bin? e vimos algumas “regras” para determinar a largura do bin

Vamos agora considerar outras técnicas para estimar densidades: KDE e GMM

Vamos exemplificar com uma amostra de elipticidades de imagens de galáxias de baixo brilho superficial na amostra do DR1 do S-PLUS (Mendes de Oliveira et al. 2019). A elipticidade é definida como \((1-B/A)\), onde \(A\) e \(B\) são o eixo maior e menor de uma galáxia

# dados:
tabela = read.table(file="ell.dat", header=TRUE)

# dimensões do arquivo:
dim(tabela)

[1] 241   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))

# histograma:
hist(ell,xlab = '1-B/A', col='salmon',main='')

# distribuições normalizadas: histograma com freq = FALSE
# a integral do histograma é igual a 1
hist(ell,xlab = '1-B/A', col='salmon',main='',freq=FALSE)

---

---

1.2 estimativa de densidades: KDE (kernel density estimation)

KDE é uma técnica que permite representar uma distribuição de forma contínua (ao contrário da representação discreta de um histograma)

KDE aplica-se a espaços de dados de qualquer dimensão \(D\)

KDE estima a densidade em um dado ponto a partir de uma média de pontos próximos usando uma função kernel

densidade (pdf) num ponto \(x\) (multi-dimensional): \[ {\hat f(x;h)} = \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

• \(N\): número de pontos

• \(D\): dimensão do espaço

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

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

---

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:

# 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))
# plota uma gaussiana
plot (density(0, bw = 1), xlab = "", main = "kernels com bw = 1")
# plota os outros kernels
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 usar kernels para examinar a distribuição de elipticidades:

as figuras abaixo têm valores diferentes de h (band width):

# efeito da largura de banda
par(mfrow = c(2,2))
h=c(1/80,1/40,1/20,bw.SJ(ell))
## a largura de banda SJ (Sheather & Jones 1991) é geralmente considerada a 
# escolha default

for(j in 1:4){
plot(density(ell, bw = h[j]), main = "mesmo h, 7 kernels diferentes")
for(i in 1: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.

nesse caso, determina-se \(h\) encontrando o valor que minimiza a expressão: \[ E\Big[ \int [\hat f_h(x) - f(x)]^2 dx \Big] \]

esta otimização é equivalente a minimizar \[ J(h) = \int \hat f_h(x)^2 dx - \frac{2}{N} \sum_{i=1}^N \log {\hat f}_{-i}(x_i;h)\] onde \({\hat f}_{-i}(x_i;h)\) é o valor da função sem o ponto \(x_i\)

# 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,xlab='h',ylab='J(h)')

Warning: `qplot()` was deprecated in ggplot2 3.4.0.

# optimal value: (mínimo)
a = optimize(J,interval=c(.01,1))
a

$minimum
[1] 0.06806629

$objective
[1] -1.493732

hopt = a$minimum
hopt

[1] 0.06806629

hopt=round(hopt,3)

# visualização:
par(mfrow = c(1,1))
y = density(ell, bw = hopt)
plot(y,xlab = '1-B/A', type = "n", main = sprintf(paste('h = ',hopt)))
polygon(y, col = "wheat")

comentários:

• a estimativa de \(h\) por validação cruzada é bastante robusta

• o uso de kernels pemite modelar as distribuições discretas como funções contínuas

• tome cuidado com soluções não físicas (como elipticidades negativas na figura acima)

• note que em alguns casos \(h\) deve ser estimado a partir de considerações físicas (por exemplo o tamanho das estruturas que se quer detectar em um certo sistema físico).

---

---

Vamos agora considerar o KDE em duas dimensões: função kde2d()

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)
suppressMessages(library(KernSmooth))
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,3))
hist(rmi,xlab='r-i',main='',col = 'salmon')
hist(gmr,xlab='g-r',main='',col = 'salmon')
plot(rmi,gmr,xlab='r-i',ylab='g-r',col = 'salmon',asp=1)

# 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')

# note a chamada "bimodalidade" na distribuição de cores

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.15)
image(f,xlab='r-i',ylab='g-r',main='h = 0.15')
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

veja como a escolha de \(h\) determina o que você vê na figura: para \(h\) pequeno a variância é muito grande, para \(h\) grande a suavização é muito forte

---

---

1.2 Gaussian Mixture Models (GMM)

nesse tipo de modelo uma distribuição (pode ser multidimensional) é modelada como uma soma de gaussianas:

• a densidade de \(n\) pontos numa distribuição modelada por \(M\) gaussianas é dada por \[\rho (x) = n \sum_{i=1}^M w_i N(x|\mu_i,\Sigma_i),\] com \(\sum_i w_i =1\) e onde \(\mu_i\) e \(\Sigma_i\) são a média e a matriz de covariância da \(i\)-ésima gaussiana

• \(w_i\) é o peso da \(i\)-ésima gaussiana

• para se obter os parâmetros das gaussianas pode-se maximizar a verossimilhança, que é proporcional a \(\rho\)

• um algoritmo muito eficiente para isso é o de Maximização de Expectativa (EM- expectation-maximization)

• o EM é um algoritmo iterativo cujo resultado depende muito da inicialização de parâmetros

vamos usar o pacote mclus para ilustrar este algoritmo com os dados de distribuição de elipticidades de galáxias; esse algoritmo implementa diversos métodos com EM, inicializados com modelos de agrupamento hierárquico; o melhor modelo é escolhido usando BIC: \(BIC = k ln(n) - 2 ln(\hat L)\)
(obs: a rotina retorna -BIC!)

a biblioteca mclus pode ser usada para estimativa de densidades, análise de agrupamento e classificação:

• estimativa de densidades: modelo <- densityMclust(dados)

• clustering: modelo <- Mclust(dados)

• classificação: mod2 <- MclustDA(dados)

vamos ilustrar com a amostra de distribuição de elipticidades:

library(mclust)

Package 'mclust' version 6.1.2
Type 'citation("mclust")' for citing this R package in publications.

set.seed(123)

# vamos usar a função densityMclust()
modelo = densityMclust(ell)

# resultado
summary(modelo)

------------------------------------------------------- 
Density estimation via Gaussian finite mixture modeling 
------------------------------------------------------- 

Mclust E (univariate, equal variance) model with 2 components: 

 log-likelihood   n df    BIC      ICL
       84.05959 241  4 146.18 76.89866

# o melhor modelo no caso tem 2 gaussianas
# acima ICL (Integrated Complete-data Likelihood) é uma outra métrica para 
# comparar modelos

summary(modelo, parameters = TRUE)

------------------------------------------------------- 
Density estimation via Gaussian finite mixture modeling 
------------------------------------------------------- 

Mclust E (univariate, equal variance) model with 2 components: 

 log-likelihood   n df    BIC      ICL
       84.05959 241  4 146.18 76.89866

Mixing probabilities:
        1         2 
0.4676829 0.5323171 

Means:
        1         2 
0.2502545 0.5200363 

Variances:
         1          2 
0.01272305 0.01272305 

# aqui temos também os parâmetros das gaussianas

# visualização:
plot(modelo, what="density",xlab="1-B/A",lwd=3)
rug(ell)

o modelo escolhido tem 2 gaussianas, mas como varia o BIC com o número de gaussianas?

# note que essa rotina retorna -BIC, de modo que quanto maior este valor, 
# mais provável o modelo

# BIC para vários modelos:
BIC <- mclustBIC(ell)

# melhores modelos:
summary(BIC)

Best BIC values:
            E,2        E,3        E,1
BIC      146.18 144.740868 143.506997
BIC diff   0.00  -1.439115  -2.672987

plot(BIC)

# o melhor modelo tem 2 gaussianas
# E e V são modelos diferentes de EM

vamos agora ilustrar GMM com a distribuição cor-cor:

# vamos colocar as cores em um único objeto:
X = cbind(rmi,gmr)

# distribuição de densidades
modelo = densityMclust(X)

summary(modelo)

------------------------------------------------------- 
Density estimation via Gaussian finite mixture modeling 
------------------------------------------------------- 

Mclust VVE (ellipsoidal, equal orientation) model with 3 components: 

 log-likelihood    n df      BIC      ICL
       1541.051 1025 15 2978.115 2424.286

# o melhor modelo tem 3 gaussianas

# parâmetros das gaussianas:
summary(modelo, parameters = TRUE)

------------------------------------------------------- 
Density estimation via Gaussian finite mixture modeling 
------------------------------------------------------- 

Mclust VVE (ellipsoidal, equal orientation) model with 3 components: 

 log-likelihood    n df      BIC      ICL
       1541.051 1025 15 2978.115 2424.286

Mixing probabilities:
        1         2         3 
0.2562417 0.3068923 0.4368659 

Means:
         [,1]      [,2]      [,3]
rmi 0.3351472 0.1192325 0.1829171
gmr 0.7021600 0.3760590 0.4566548

Variances:
[,,1]
            rmi         gmr
rmi 0.003849905 0.002757176
gmr 0.002757176 0.007369840
[,,2]
            rmi         gmr
rmi 0.016903791 0.008049297
gmr 0.008049297 0.027179883
[,,3]
            rmi         gmr
rmi 0.004876661 0.005890428
gmr 0.005890428 0.012396645

# visualização:
plot(modelo, what="density",xlab='(r-i)',ylab='(g-r)',lwd=2,col='blue')

# vamos ver o BIC para vários modelos:
BIC <- mclustBIC(X)
plot(BIC)

# a figura ilustra resultados com vários modelos de inicialização

# melhores modelos:
summary(BIC)

Best BIC values:
            VVE,3       VVV,3      VEE,3
BIC      2978.115 2976.686030 2966.05180
BIC diff    0.000   -1.428725  -12.06295

# o melhor modelo tem 3 gaussianas

preste atenção em como a forma da distribuição depende do método adotado para modelar os dados!

---

---

2. análise de agrupamento

objetivo: encontrar grupos (clusters) no espaço de dados

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,…

---

---

para identificar grupos, considera-se as distâncias entre os dados

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

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

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

• mahalanobis: \(d_{ij} = \Big[(\mathbf{x}_i - \mathbf{x}_j)^T.\mathbf{C}^{-1}.(\mathbf{x}_j-\mathbf{x}_j) \Big]^{1/2}\), onde \(C\) é a matriz de covariância

---

---

2.1 algoritmo K-means:

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

• custo: 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’)

Note que:

• K deve ser dado a priori

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

• os resultados podem depender do tipo de distância escolhida

Vamos considerar inicialmente, para ilustrar, dados simulados em duas dimensões:

# reprodutibilidade
set.seed(123)

# vamos gerar coordenadas (x,y) para dois conjuntos de dados
# 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))

# vamos juntar os dois conjuntos de dados em um único conjunto:
dados = rbind(grupo1, grupo2)

#visualização
par(mfrow = c(1,1))
plot(dados, pch = 16, xlab = "x", ylab = "y", xlim=c(-5,5),ylim=c(-5,5))

# a função kmeans do R permite rodar o algoritmo com várias inicializações 
# aleatórias diferentes, retornando o resultado de menor custo J
# 2 grupos e 20 reinicializações:
kmres = kmeans(dados,2, nstart = 20)

# centros dos grupos: compare com mu1 e mu2
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="yellow", 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", 'yellow'), 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, que mistura os grupos:
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)

plot(dados, pch = 16, xlab = "x", ylab = "y", xlim=c(-5,5),ylim=c(-5,5))

# 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 agora aplicar o K-means ao diagrama cor-cor de uma amostra de galáxias do S-PLUS

Para ilustrar os dados vamos usar a rotina kernel2d, que ajusta um kernel em 2 dimensões, mais a função image:

# kernel2d
par(mfrow = c(1,1))
library(MASS)
dados_sm = kde2d(rmi, gmr, h=0.075,n=500)

# visualização do resultado com image()
image(dados_sm, col=grey(13:0/15), xlab='r-i',ylab='g-r',main='galáxias do S-PLUS')

pré-processamento
em problemas de aprendizado de máquina a escala dos dados pode ser importante e, em muitos casos, é recomendável se fazer algum tipo de pré-processamento dos dados, para que todas as variáveis fiquem mais ou menos com a mesma escala.

os dois tipos de pré-processamento mais comuns são a) renormalizar cada variável para que fique entre 0 e 1 ou entre -1 e 1 ou, b) subtrair a média e dividir pelo desvio padrão da variável (“padronização”)

# dados
dados=cbind(rmi,gmr)

# padronização das variáveis
# o default do scale é: scale(x, center = TRUE, scale = TRUE):
# subtrai a média e divide pelo desvio padrão
dados_scale = scale(dados)

# como o custo varia com o número de grupos?
# vamos considerar 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)

2.2 análise com GMM

vamos agora modelar a distribuição de cores com o algoritmo GMM:

# vamos colocar as cores em um único objeto:
X = cbind(rmi,gmr)

# vamos ver o BIC para vários modelos:
BIC <- mclustBIC(X)
plot(BIC)

# a figura ilustra resultados com vários modelos de inicialização

# melhores modelos:
summary(BIC)

Best BIC values:
            VVE,3       VVV,3      VEE,3
BIC      2978.115 2976.686030 2966.05180
BIC diff    0.000   -1.428725  -12.06295

# o melhor modelo tem 3 gaussianas
# NOTE QUE ESTES MODELOS SÃO EXATAMENTE OS MESMOS OBTIDOS NA ESTIMATIVA DE 
# DENSIDADES ACIMA!

mod1 <- Mclust(X, x = BIC)
summary(mod1, parameters = TRUE)

---------------------------------------------------- 
Gaussian finite mixture model fitted by EM algorithm 
---------------------------------------------------- 

Mclust VVE (ellipsoidal, equal orientation) model with 3 components: 

 log-likelihood    n df      BIC      ICL
       1541.051 1025 15 2978.115 2424.286

Clustering table:
  1   2   3 
276 232 517 

Mixing probabilities:
        1         2         3 
0.2562417 0.3068923 0.4368659 

Means:
         [,1]      [,2]      [,3]
rmi 0.3351472 0.1192325 0.1829171
gmr 0.7021600 0.3760590 0.4566548

Variances:
[,,1]
            rmi         gmr
rmi 0.003849905 0.002757176
gmr 0.002757176 0.007369840
[,,2]
            rmi         gmr
rmi 0.016903791 0.008049297
gmr 0.008049297 0.027179883
[,,3]
            rmi         gmr
rmi 0.004876661 0.005890428
gmr 0.005890428 0.012396645

plot(mod1, what = "classification")

note como a classe em vermelho se sobrepõe com a verde!
GMM não é uma boa técnica de clustering para esses dados!

---

---

3. PCA - Análise de Componentes Principais

Vocês precisam conhecer PCA!

Se temos dados com N atributos (dimensões), será que conseguimos representá-los com menos variáveis?

Para que serve?

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

• visualização de dados multidimensionais

como funciona? PCA cria um novo sistema de coordenadas ortonormal no espaço de dados:

• 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\)

• 2. 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\)

• 3. no sub-espaço perpendicular a \(PC_1\) e \(PC_2\) identifique…

• 4. repita o procedimento até \(PC_D\), onde \(D\) é a dimensão do espaço de dados

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

isso significa que os dados podem ser aproximadamente descritos com um número \(d << D\) de componentes.

---

---

matematicamente:

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

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

um objeto é um ponto no espaço de dados \(D\)-dimensional

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} ~~~~~~~~~j=1,...,D \] 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\): \[ C\mathbf{a}_j = \lambda_j \mathbf{a}_j\]

o autovetor de \(C\) com o j-ésimo maior autovalor \(\lambda_j\) é 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 \(\lambda_j\) 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

---

---

Exemplo: vamos analisar uma amostra contendo 100 espectros de galáxias

Esses dados vêm da síntese espectral de galáxias observadas pelo SDSS feita com o código starlight por Werle et al. (2019, MNRAS, 483, 2382)

Entre outras coisas, o starlight produz um modelo para o espectro de uma galáxia

Aqui selecionei 100 galáxias ao acaso, e vamos analisar os modelos para o contínuo e linhas de absorção dos espectros, entre 3500 e 5999\(\overset{\circ}{A}\) (2500 valores de comprimento de onda)

Os espectros estão normalizados em 5500\(\overset{\circ}{A}\)

Supõe-se que o contínuo e as linhas de absorção sejam formadas na fotosfera das estrelas, enquanto as linhas de emissão, não incluidas nesses dados, são produzidas em regiões gasosas ionizadas

Note que, nesse exemplo, um espectro é um ponto em um espaço de 2500 dimensões! E vamos mostrar que podem ser bem descritos com muito poucas dimensões!!!

leitura dos espectros:

# nessa tabela cada linha corresponde a um comprimento de onda e cada coluna 
# ao fluxo nesse comprimento de onda de uma galáxia diferente
espectros = read.table('spec100c.dat',header=T, fill=T)
dim(espectros)

[1] 2500  101

#2500  101

#comprimentos de onda: a primeira coluna
lambdaf = espectros[,1]

#fluxos: todas as colunas, menos a primeira
fluxos = espectros[,-1]

dim(fluxos)

[1] 2500  100

# 100 objetos, 2500 fluxos por objeto

Exemplos de alguns espectros:

par(mfrow = c(2,3))
# exemplo de espectros 
plot(lambdaf,fluxos[,1],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')
plot(lambdaf,fluxos[,7],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')
plot(lambdaf,fluxos[,12],
type='l',xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')
plot(lambdaf,fluxos[,28],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')
plot(lambdaf,fluxos[,47],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')
plot(lambdaf,fluxos[,65],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon')

todos os espectros:

par(mfrow = c(1,1))
plot(lambdaf,fluxos[,1],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='fluxo',col='salmon',ylim=c(0,2))
for(i in 2:ncol(fluxos)){
lines(lambdaf,fluxos[,i],col='salmon')
}

Vamos agora calcular o PCA dos fluxos. Para isso vamos considerar os dados representados como uma matriz onde cada linha é uma galáxia e as colunas são os fluxos nos 2500 comprimentos de onda: essa é a matriz transposta da matrix de fluxos acima

dados_pca = t(fluxos)

A coluna correspondente a 5500\(\overset{\circ}{A}\) tem o valor 1 para todas as galáxias, devido à normalização dos espectros. O PCA não gosta disso e vamos, então, remover essa coluna da análise e do vetor com os comprimentos de onda, lambdaf:

#para ver qual é ela:
inorm = which(lambdaf == 5500)
inorm

[1] 2001

# removendo ela:
dados_pca = dados_pca[,-inorm]

dim(dados_pca)

[1]  100 2499

# vamos remover esse valor de lambdaf também:
lambdaf=lambdaf[-inorm]

Para o PCA vamos escalonar as variáveis: média zero e variância unitária

PCA em R é feito com a função prcomp()

pca.fluxos = prcomp(dados_pca, retx=TRUE, center=TRUE, scale=TRUE)

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

# autovalores
lambda = sd^2

# fraçao da variância explicada
ve = round(cumsum(lambda)/sum(lambda),digits = 3)
ve[1:5]

[1] 0.892 0.937 0.966 0.976 0.984

# vemos que a primeira componente sozinha explica 89% da variância e as duas 
# primeiras 94%!

# visualização: scree plot - variância associada a cada componente
screeplot(pca.fluxos, main=" scree plot")

# isso significa que 94% da variância do espaço espectral de 2499 dimensões 
# está contida em um plano definido pelas duas primeiras componentes principais!

# autovetores
autovec = pca.fluxos$rotation

# componentes principais
pc = pca.fluxos$x

# 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')

# PCA: visualização do espaço de dados a partir de um referencial especial:
par(mfrow = c(2,2))
plot(pc[,1], pc[,2], xlab="PC 1", ylab="PC 2",col='blue')
plot(pc[,1], pc[,3], xlab="PC 1", ylab="PC 3",col='blue')
plot(pc[,2], pc[,3], xlab="PC 2", ylab="PC 3",col='blue')
plot(pc[,3], pc[,4], xlab="PC 3", ylab="PC 4",col='blue')

# note que estes gráficos não estão em escala: compare os intervalos de 
# valores nos dois eixos

reconstrução dos espectros:

#https://gist.github.com/menugget/7689201
#This function reconstructs a data set using a defined set of principal components.
#arguments "pca" is the pca object from prcomp, "pcs" is a vector of 
# principal components
#to be used for reconstruction (default includes all pcs)
prcomp.recon <- function(pca, ncomp){
  pcs <- seq(1:ncomp)
  recon <- as.matrix(pca$x[,pcs]) %*% t(as.matrix(pca$rotation[,pcs]))
  if(pca$scale[1] != FALSE){
    recon <- scale(recon , center=FALSE, scale=1/pca$scale)
  }
  if(pca$center[1] != FALSE){
    recon <- scale(recon , center=-pca$center, scale=FALSE)
  }
  recon
}

# vamos ilustrar a reconstrução do espectro 12 usando de 1 a 4 componentes>
par(mfcol=c(2,2))
obj = 12
for(nComp in 1:4){
Xhat = prcomp.recon(pca.fluxos,nComp)
plot(lambdaf,Xhat[obj,],type='l',
xlab=expression(paste(lambda,' (',ring(A),')')),ylab='flux',col='blue',
ylim=c(0.5,1.7))
lines(lambdaf,dados_pca[obj,],col='red')
legend('topright',legend = paste(nComp),bty='n')
}

---

---

4. redução de dimensionalidade

O objetivo da redução de dimensionalidade é reduzir o número de variáveis, ou a dimensão, de um conjunto de dados, procurando manter a estrutura intrínseca dos dados.

Dito de outro modo, queremos obter uma representação dos dados num espaço de dimensão menor, retendo o máximo de informações.

Muitas vezes há redundância nas variáveis, e a redução de dimensionalidade então procura representar essencialmente a mesma informação com menos variáveis

Um ponto importante é que os procedimentos de redução de dimensionalidade criam novas variáveis que procuram preservar ao máximo a estrutura dos dados; mas nem sempre as novas variáveis são interpretáveis

Muitas vezes o espaço definido por essas novas variáveis é chamado de espaço latente.

Para ilustrar alguns algoritmos de redução de dimensionalidade vamos considerar uma amostra extraída do S-PLUS/DR4 com classificação dos objetos em 3 classes:
0: QSO, 1:estrelas, 2: galáxias

Nessa amostra todos os objetos têm fotometria em todas as 12 bandas do S-PLUS e magnitudes 15 < r_auto < 21.

Note que alguns objetos dessa amostra podem ter sido mal-classificados pois não foi feita nenhuma seleção com base na probabilidade de um objeto ser de uma certa classe.

A tabela contém 3000 objetos, 1000 em cada uma das 3 classes.

# leitura da amostra
dados = read.table(file='qso_est_gal.dat',head = TRUE)
colnames(dados)

 [1] "uJAVA"  "F378"   "F395"   "F410"   "F430"   "g_auto" "F515"   "r_auto"
 [9] "F660"   "i_auto" "F861"   "z_auto" "class" 

dim(dados)

[1] 3000   13

# vamos definir uma variável, *class*, que vai conter os tipos dos objetos:
class  = dados[,13]

head(class)

[1] 1 1 0 2 1 1

# vamos fazer uma função para plotar dados, com "transparência" nos pontos 
# (controlada pelo parâmetro alpha: alpha = 1- opaco, alpha = 0 - totalmente 
# transparente)
plota_transp = function(x,y,nome,nomex,nomey){
plot(x[class == 0],y[class == 0],pch=20,col=alpha('black', 0.2),main=nome,
xlim = c(min(x),max(x)),asp = 1,xlab=nomex,ylab=nomey,cex.lab=1.5)
legend("topleft",legend = levels(factor(class)),pch = 19,
col = factor(levels(factor(class))))
points(x[class == 1], y[class == 1],col=alpha('red', 0.2))
points(x[class == 2], y[class == 2],col=alpha('green', 0.2))
}

# vamos ilustrar esta função com um diagrama cor-cor
par(mfrow = c(1,1))
plota_transp(dados$g_auto-dados$r_auto,dados$i_auto-dados$z_auto,
'diagrama cor-cor','(g-r)','(i-z)')

Esta amostra tem dimensão 12: as magnitudes do S-PLUS. Para verificar se há correlação entre as variáveis podemos calcular a matriz de correlação da amostra:

library(corrplot)

corrplot 0.95 loaded

library(RColorBrewer)
M <-cor(dados[,1:12],method = 'spearman')
par(mfrow = c(1,1))
corrplot(M, type="upper", order="hclust",col=brewer.pal(n=8, name="RdYlBu"))

Examinando-se o gráfico, pode-se verificar que, de fato, há variáveis bastante correlacionadas, o que sugere que um espaço latente com \(n << 12\) pode ser uma boa representação dos dados dessa amostra.

---

---

4.1 PCA

O primeiro método que consideraremos é o PCA, que é um método linear:

# vamos escalonar as variáveis: média zero e variância unitária
pca.fot = prcomp(dados[,1:12], retx=TRUE, center=TRUE, scale=TRUE)

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

# autovalores
lambda = sd^2

# fraçao da variância explicada
ve = round(cumsum(lambda)/sum(lambda),digits = 3)
ve

 [1] 0.827 0.940 0.957 0.970 0.980 0.990 0.996 0.998 0.999 0.999 1.000 1.000

# 0.827 0.940 0.957 0.970 0.980 0.990 0.996 0.998 0.999 0.999 1.000 1.000
# vemos que as duas primeiras componentes explicam 94% da variância

# visualização: scree plot - variância associada a cada componente
screeplot(pca.fot, main=" scree plot")

# componentes principais
pc = pca.fot$x

# visualização: duas primeiras componentes
par(mfrow = c(1,1))
plota_transp(pc[,1], pc[,2],'','PC1','PC2')

par(mfrow = c(2,2))
plota_transp(pc[,1], pc[,2],'','PC1','PC2')
plota_transp(pc[,1], pc[,3],'','PC1','PC3')
plota_transp(pc[,2], pc[,3],'','PC2','PC3')
plota_transp(pc[,3], pc[,4],'','PC3','PC4')

O PCA acima foi rodado escalonando as variáveis para ficarem com média nula e variância unitária.

---

---

4.2 Uniform manifold approximation and projection (UMAP)

UMAP produz uma representação de baixa dimensão de dados de alta dimensão que preserva a estrutura relevante, mantendo uma estrutura topológica a mais próxima possível dos dados originais de alta dimensão.

Vamos considerar projeções em 2 e 3 dimensões:

library(umap)

# UMAP em 2 dimensões
umap_fot = umap(dados, n_components = 2, random_state = 15) 

layout <- umap_fot[["layout"]] 

layout <- data.frame(layout) 

plot(layout,col=factor(class),xlab='X1',ylab='X2', pch=16, cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(class)),
       pch = 19,
       col = factor(levels(factor(class))))

# UMAP em 3 dimensões
umap_fot = umap(dados, n_components = 3, random_state = 15) 

layout <- umap_fot[["layout"]] 

layout <- data.frame(layout) 

plot(layout[,1],layout[,2],col=factor(class),xlab='X1',ylab='X2', pch=16, 
cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(class)),
       pch = 19,
       col = factor(levels(factor(class))))

plot(layout[,1],layout[,3],col=factor(class),xlab='X1',ylab='X3', pch=16, 
cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(class)),
       pch = 19,
       col = factor(levels(factor(class))))

plot(layout[,2],layout[,3],col=factor(class),xlab='X2',ylab='X3', pch=16, 
cex=2, cex.lab=1.5) 
legend("topleft",
       legend = levels(factor(class)),
       pch = 19,
       col = factor(levels(factor(class))))

Problemas de redução de dimensionalidade envolvem cálculos de distâncias no espaço dos dados. Como variáveis diferentes podem ter um intervalo de variação muito diferente, pode ser importante re-escalonar os dados antes da análise.

O espaço reduzido é chamado espaço latente.

Note que a morfologia na distribuição dos pontos no espaço latente obtida com diferentes métodos é diferente e pode ser explorada para descobrir relações ocultas entre os dados ou para se fazer classificação de objetos.

---

---

5. Exercícios

Inicialize a semente de números aleatórios com seu número USP.

1. O arquivo spec100cl.dat contém modelos no referencial de repouso para as mesmas 100 galáxias do arquivo spec100c.dat, mas incluindo tanto o contínuo e linhas de absorção, quanto linhas de emissão.

• 1.1. Examine os dados fazendo figuras com alguns espectros.

• 1.2. Verifique quais são as linhas de emissão mais importantes na região coberta por esses espectros.

• 1.3. Faça uma análise desses dados com PCA e verifique como muda a variância associada às primeiras componentes principais em relação ao caso sem as linha.

• 1.4. Examine a reconstrução dos espectros nesse caso.

2. Aplique UMAP a um conjunto de dados constituido pelas 5 primeiras componentes principais determinadas no exercício anterior. Veja como o resultado varia com o parâmetro n_neighbors (considere valores 3, 5, 10 e 20). Comente sobre como a aparência dos plots de UMAP1 x UMAP2 muda.

---

---

`