aula9

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 compressão de informação

3.1 PCA: análise de componentes principais

3.2 o pacote dimRed

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

# tamanho do bin
n_silverman = ldata/bin_silverman+1
round(n_silverman)

## [1] 13

# regra de Knuth: teste com até 100 bins

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

N = length(data)

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

logp = rep(0,N)

for(M in 2:maxM){

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

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

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

a = optBINS(ell,100)

Mknuth = which(a == max(a))
Mknuth

## [1] 8

# 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

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

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 2:length(kernels)){lines(density(ell, bw = h[j], kernel = kernels[i]), col = i)}}

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

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

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

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

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

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

## $minimum
## [1] 0.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:

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 5.4.10
## 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 é 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:

# 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. Compressão de informação

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

4.1 PCA - Análise de Componentes Principais

vocês precisam conhecer PCA!

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

# gráficos em densidade
par(mfrow = c(2,2))
smoothScatter(pc[,1],pc[,2],nrpoints=0,add=FALSE, xlab='pc1', ylab='pc2', main = '',asp = 1)
smoothScatter(pc[,1],pc[,3],nrpoints=0,add=FALSE, xlab='pc1', ylab='pc3', main = '',asp = 1)
smoothScatter(pc[,2],pc[,3],nrpoints=0,add=FALSE, xlab='pc2', ylab='pc3', main = '',asp = 1)
smoothScatter(pc[,3],pc[,4],nrpoints=0,add=FALSE, xlab='pc3', ylab='pc4', main = '',asp = 1)

>

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.2 o pacote dimRed

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

suppressMessages(library(dimRed))

para saber quais são os métodos disponíveis:

dimRedMethodList()

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

o manual do pacote contém os hiperparâmetros necessários para rodar cada método

vamos preparar os dados colocando um objeto por linha (ao invés de coluna):

dados <- t(as.matrix(fluxos))
# normalização
#dados <- scale(dados, center = TRUE, scale = TRUE)
dim(dados)

## [1]  100 2500

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

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

vamos fazer um teste com LLE (locally linear embbeding):

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

• vizinhança: k-vizinhos mais próximos

• algoritmo não-supervisionado

• o algoritmo é iterativo:

• 1. para cada ponto determina-se um conjunto de pesos W que melhor o reconstroi a partir de seus k-vizinhos mais próximos (ajuste de um hiperplano)- esses pesos codificam a geometria local

• 2. com os pesos fixos, determina-se um conjunto de dados de menor dimensão que preserva as relações de vizinhança desses pesos; os pesos W são obtidos via álgebra linear

• 3. a análise depende do número de vizinhos de cada ponto: se o número de vizinhos é pequeno, o mapeamento não vai refletir nenhuma propriedade global, enquanto que se for muito grande vai perder seu caráter não-linear e se comportar como PCA

# este algoritmo depende de um hiperparâmetro: o número de vizinhos
# podemos encontrar um valor 'ótimo' com um código do tipo:
kk <- floor(seq(1, 15, 1))
# embedding:
emb <- lapply(kk, function(x) embed(dados, "LLE", knn = x,.mute='output'))
## Qualidade dos embeddings
qual <- sapply(emb, function(x) quality(x, "Q_local"))

## Warning in rankmatrix(X, input = "dist", use): 0 outside of diagonal in distance
## matrix

## Warning in rankmatrix(X, input = "dist", use): 0 outside of diagonal in distance
## matrix

## Warning in rankmatrix(X, input = "dist", use): 0 outside of diagonal in distance
## matrix

## melhor valor para knn
ind_max <- which.max(qual)
k_max <- kk[ind_max]
k_max

## [1] 8

# solução com o melhor valor de knn
emb_LLE <- embed(dados, "LLE", knn = k_max)

## finding neighbours
## calculating weights
## computing coordinates

plot(emb_LLE, type = "2vars")

vamos agora testar outros métodos desse pacote: PCA, kernel PCA, e Isomap

# PCA
emb_PCA <- embed(dados, "PCA")
plot(emb_PCA, type = "2vars")

# kPCA  kernel PCA
# extensão não-linear do PCA usando kernels
emb_kPCA <- embed(dados, "kPCA")
plot(emb_kPCA, type = "2vars")

# Isomap
# usa a distância geodésica e não a euclidiana para promover o embedding
# o método tem dois hiperparâmetros: knn e ndim, o número de dimensões do embedding (default = 2)
# vamos encontrar um valor 'ótimo' para knn:
kk <- floor(seq(10, 60, 1))
# embedding:
 emb <- lapply(kk, function(x) embed(dados, "Isomap", knn = x, .mute = c('message','output')))
## Qualidade dos embeddings
qual <- sapply(emb, function(x) quality(x, "Q_local"))
## melhor valor para knn
ind_max <- which.max(qual)
k_max <- kk[ind_max]
k_max

## [1] 38

# solução com o melhor valor de knn
emb_Isomap <- embed(dados, "Isomap", knn = k_max)
plot(emb_Isomap, type = "2vars")

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

do manual do dimRed:

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

vamos avaliar a qualidade da redução de dimensionalidade dos vários métodos:

# Q_local
quality(emb_PCA,'Q_local')

## [1] 0.69815

quality(emb_kPCA,'Q_local')

## [1] 0.528593

quality(emb_LLE,'Q_local')

## [1] 0.5938881

quality(emb_Isomap,'Q_local')

## [1] 0.7536438

# ###

# Q_global
quality(emb_PCA,'Q_global')

##       14 
## 0.408004

quality(emb_kPCA,'Q_global')

##        25 
## 0.1927997

quality(emb_LLE,'Q_global')

##        18 
## 0.2772082

quality(emb_Isomap,'Q_global')

##        12 
## 0.4205566

# ###

# R_NX
quality(emb_PCA,'mean_R_NX')

## [1] 0.9505007

quality(emb_kPCA,'mean_R_NX')

## [1] 0.4311157

quality(emb_LLE,'mean_R_NX')

## [1] 0.6432584

quality(emb_Isomap,'mean_R_NX')

## [1] 0.9641862

neste exemplo o Isomap conseguiu um nível de compressão maior que PCA

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

---

---

Exercícios

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, KDE e GMM.

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

• 2.5. Use o pacote dimRed e compare a compressão de dados obtida com PCA e kPCA.

---

---

`