aula de hoje: análise de dados não paramétrica

1 estimativa de densidades:

1.1 histogramas

1.2 KDE kernel density estimation

1.3 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 ICA: Independent component analysis

4.3 t-SNE: t-distributed stochastic neighbor embedding

4.4 UMAP: Uniform manifold approximation and projection

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

1.1 estimativa de densidades: histogramas

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

considere o problema: temos um conjunto de N dados que queremos representar com um histograma:
qual é o tamanho “ideal” do bin?

existe um conjunto de regras para determinar a largura do bin
seja \(\sigma_s\) a dispersão da amostra e \(\Delta\) a largura do bin; então:

  • regra de Scott: \[ \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 assumindo uma verossimilhança multinomial e um prior não-informativo
    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 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))
# vamos usar nesta figura 40 bins
hist(ell,xlab = '1-B/A', col='salmon',main='',breaks=40)

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

número “ótimo” de bins de acordo com as várias regras:

# regra de Scott:
sigma = sd(ell)
# tamanho do bin
bin_scott = 3.5*sigma/length(ell)^(1/3)
bin_scott
## [1] 0.09897846
# 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
# versão robusta do desvio padrão de uma gaussiana
sig_G = 0.7413*(quantile(ell,0.75,names = FALSE) - quantile(ell,0.25,names = FALSE))
# tamanho do bin
bin_fd = 3.5*sig_G/length(ell)^(1/3)
bin_fd
## [1] 0.115858
# número de bins
n_fd = ldata/bin_fd+1
round(n_fd)
## [1] 7
# regra de Silverman:
# tamanho do bin
bin_silverman = 1.06*sigma/length(ell)^(1/5)
bin_silverman
## [1] 0.06228461
# número de bins
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
}
return(logp)
}

a = optBINS(ell,100)

Mknuth = which(a == max(a))
Mknuth
## [1] 8
# distribuição de probabilidades de M
par(mfrow = c(1,1))
plot(seq(2,100,1),a[2:100],type='l',xlab='M', ylab='log(prob)',col='red',lwd=3)
abline(v=8,lwd=2)

# número de bins dos vários métodos:
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=0.9,bty='n')
plot(cut(ell, breaks = n_fd),col='green')
legend('topleft', 'Freedman-Diaconis',cex=0.9,bty='n')
plot(cut(ell, breaks = n_silverman),col='darkgoldenrod')
legend('topleft', 'Silverman',cex=0.9,bty='n')
plot(cut(ell, breaks = Mknuth),col='blue')
legend('topleft', 'Knuth',cex=0.9,bty='n')

note que a escolha do tamanho dos bins não é óbvia pois, a rigor, o resultado deve depender da forma intrínseca da distribuição, que não se conhece!

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:

# 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.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# optimal value:
a = optimize(J,interval=c(.01,1))
a
## $minimum
## [1] 0.06811105
## 
## $objective
## [1] -1.493788
hopt = a$minimum
hopt
## [1] 0.06811105
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')

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.0.1
## 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
    1. inicialize, escolhendo K objetos ao acaso como centro
    1. 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:

# 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\)
    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\), 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.92 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.

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, é importante re-escalonar os dados antes da análise.

Vamos então criar novas variáveis, escalonadas entre 0 e 1:

# reescalonamento usando a função scale():
maxs <- apply(dados[,1:12], 2, max) 
mins <- apply(dados[,1:12], 2, min)
dados_scaled <- as.data.frame(scale(dados[,1:12], center = mins, scale = maxs - mins))

dados = dados_scaled

4.2 Independent component analysis (ICA)

ICA decompõe um sinal multivariado em sinais não gaussianos estatisticamente independentes.

Uma aplicação típica é o “cocktail party problem”: temos várias pessoas conversando ao mesmo tempo em uma sala e o algoritmo separa as conversas.

library(ica)
ica_fot <- icafast(dados[,1:12], 2, center = TRUE, maxit = 100, tol = 1e-6)

par(mfrow = c(1,1))
ICA1 = ica_fot$Y[, 1]
ICA2 = ica_fot$Y[, 2]
plota_transp(ICA1, ICA2,'','ICA1','ICA2')

4.3 t-distributed stochastic neighbor embedding (t-SNE)

t-SNE é uma técnica de redução de redução de dimensionalidade não-linear útil para visualização dos dados em 2 ou 3 dimensões.

Ela modela cada ponto em altas dimensões como um ponto em 2 ou 3 dimensões, tal que objetos similares são modelados como pontos próximos e objetos diferentes como pontos distantes, com alta probabilidade.

t0 = Sys.time()
library(Rtsne)
set.seed(123)
tsne_fot <- Rtsne(dados[,1:12],pca = FALSE, perplexity = 10, theta = 0.0)
TSNE1 = tsne_fot$Y[, 1]
TSNE2 = tsne_fot$Y[, 2]

par(mfrow = c(1,1))
plota_transp(TSNE1, TSNE2,'','TSNE1','TSNE2')

Sys.time()-t0
## Time difference of 2.342288 mins

O t-SNE é um método estocástico e cada vez que se roda se obtém distribuições diferentes no espaço latente, de modo que, para assegurar a reprodutibilidade, é importante se definir uma semente para inicializar o gerador de números aleatórios do algoritmo.

Outro ponto importante é que muitos algoritmos para redução de dimensionalidade dependem de certos parâmetros que precisam ser escolhidos.

O t-SNE, por exemplo, tem um parâmetro, ‘perplexity’, que controla o balanço entre propriedades locais e globais dos dados: por exemplo, se a ‘perplexity’ é baixa, muitos grupos pequenos se formam. Vamos ilustrar como fica os resultados desse algoritmo para vários valores desse parâmetro:

t0 = Sys.time()
perplexity <- seq(2, 26, 8)
set.seed(123)
par(mfrow = c(2,2))
for(i in 1:length(perplexity)){
tsne_fot <- Rtsne(dados[,1:12],pca = FALSE, perplexity = perplexity[i], theta = 0.0)
TSNE1 = tsne_fot$Y[, 1]
TSNE2 = tsne_fot$Y[, 2]
plota_transp(TSNE1, TSNE2,perplexity[i],'TSNE1','TSNE2')
}

Sys.time()-t0
## Time difference of 9.188944 mins

4.4 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.

UMAP é semelhante ao t-SNE, mas preserva a estrutura mais global nos dados.

library(uwot)
## Loading required package: Matrix
set.seed(123)
umap_fot <- umap(dados[,1:12], n_neighbors = 15, min_dist = 1, spread = 5)
UMAP1 = umap_fot[, 1]
UMAP2 = umap_fot[, 2]

par(mfrow = c(1,1))
plota_transp(UMAP1, UMAP2,'Uniform manifold approximation and projection (UMAP)','UMAP1','UMAP2')

Note que a morfologia na distribuição dos pontos em um espaço de duas dimensões obtida com cada método é diferente e pode ser explorada para descobrir relações ocultas entre os dados ou para se fazer classificação de objetos.

Exercícios

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

  1. O arquivo zspec.dat contém redshifts espectroscópicos entre 0 e 0.5. Quantas estruturas você acha que tem ao longo da linha de visada? Use técnicas de modelagem da distribuição de probabilidades para abordar o problema. Considere a regra de Knuth e GMM.
  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.
  • 2.1. Examine os dados fazendo figuras com alguns espectros.
  • 2.2. Verifique quais são as linhas de emissão mais importantes na região coberta por esses espectros.
  • 2.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.
  • 2.4. Examine a reconstrução dos espectros nesse caso.
  1. 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.

`