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 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}}$ onde $q_{25}$ e $q_{75}$ são quartis 25% e 75% da distribuição

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
# tamanho do bin
bin_fd = 2*(quantile(ell,0.75,names = FALSE) - quantile(ell,0.25,names = FALSE))/length(ell)^(1/3)
bin_fd

## [1] 0.08930878

# número de bins
n_fd = ldata/bin_fd+1
round(n_fd)

## [1] 9

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

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

defina o número de grupos, K

inicialize, escolhendo K objetos ao acaso como centro

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:

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$

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$

no sub-espaço perpendicular a $PC_1$ e $PC_2$ identifique…

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.

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

Exercícios

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

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.

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.

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.

aula8

lsj

May 9, 2025

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 UMAP: Uniform manifold approximation and projection

Exercícios

1. estimativa de densidades

1.1 estimativa de densidades: histogramas

**1.2 estimativa de densidades: KDE (kernel density estimation)**

1.2 Gaussian Mixture Models (GMM)

2. análise de agrupamento

2.1 algoritmo K-means:

2.2 análise com GMM

3. PCA - Análise de Componentes Principais

4. redução de dimensionalidade

4.1 PCA

4.2 Uniform manifold approximation and projection (UMAP)

Exercícios