M todo Kernel: Estima o de Densidades e Classifica o de ...fatc/AM/kernel.pdf · Kernel (1) NÃO é o núcleo celular NÃO é o núcleo de um sistema operacional NÃO é o núcleo/espaço

Introdução Motivação Estimação de Densidades Análise Discriminante Experimentos Conclusões Referências

Método Kernel: Estimação de Densidades e Classificação dePadrões

Marcelo Rodrigo Portela Ferreira

Departamento de Estatística, UFPB

Centro de Informática, UFPE

15 de abril de 2009

Ferreira, M. R. P. DE - UFPB / CIn - UFPE

Método Kernel: Estimação de Densidades e Classificação de Padrões


Estrutura da Apresentação

1 Motivação2 Estimação de Densidades pelo Método Kernel

(i) Caso Univariado(ii) Caso Multivariado

3 Análise Discriminante Kernel4 Experimentos5 Conclusões




Kernel (1)

NÃO é o núcleo celular

NÃO é o núcleo de um sistema operacional

NÃO é o núcleo/espaço nulo de uma matriz A : {x˜

: Ax˜

= 0˜}




Kernel (2)

Uma função K : Rp → R tal que

K(x˜) ≥ 0

∫

Rp

K(x˜)dx˜

= 1

K é simétrica em torno de 0˜




Motivação

Na metodologia clássica, faz-se alguma suposição sobre a forma funcionalparamétrica dos dados

Com uma forma paramétrica imposta, tudo que resta é estimar osparâmetros através dos dados (Máxima verossimilhança, por exemplo)

Muitas vezes, a suposição acerca da forma funcional paramétrica pode sermuito restritiva ou, em alguns casos, inadequada

Abordagens não-paramétricas permitem-nos lidar com um número maiorde situações




Estimação de Densidades pelo Método Kernel

A função densidade de probabilidade é um conceito fundamental emestatística

Uma variável aleatória X com função de distribuição F é ditaabsolutamente contínua se existir uma função não negativa f tal queF (x) = P (X ≤ x) =

∫ x

−∞ f(t)dt, ∀x ∈ R. Neste caso, dizemos que f é a

função densidade de probabilidade de X e deve satisfazer∫ +∞−∞ f(x)dx = 1

Especificar a função densidade de X nos fornece uma descrição natural dasua distribuição, e permite que probabilidades associadas a X possam serencontradas através da relação

P (a < X < b) =

∫ b

a

f(x)dx, para todo a < b.




Nosso foco: estimação de funções densidade através do método kernel

Outros tipos de estimadores não-paramétricos de funções densidadeincluem: histogramas, polígonos de frequência, splines, estimadoresbaseados em séries ortogonais e estimadores baseados em verossimilhançapenalizada (Silverman 1986; Scott 1992; Simonoff 1996)

O estimador kernel pode ser pensado como uma generalização dohistograma




Histogramas

Método não-paramétrico mais antigo de estimação de densidades

Dada uma origem x0 e um comprimento de intervalo h, definimos osretângulos do histograma como sendo os intervalos[x0 + (r − 1)h, x0 + rh) para valores inteiros positivos e negativos de r

Empiricamente a idéia e contar o número de observações que estãocontidas em cada intervalo

Sem perda de generalidade, seja o intervalo [−h/2, h/2). A probabilidadede uma observação pertencer ao intervalo [−h/2, h/2) é dada porP (X ∈ [−h/2, h/2)) =

∫ h/2

−h/2f(x)dx, onde f é a densidade de X

Uma aproximação natural para a probabilidade acima éP (X ∈ [−h/2, h/2)) ≈ 1

n#{Xi ∈ [−h/2, h/2)}

Dessa forma, uma estimativa para f seria

f(x) =1

nh#{Xi ∈ [−h/2, h/2)}, ∀x ∈ [−h/2, h/2)




Este estimador não é contínuo e depende fortemente da escolha de h,conhecido como parâmetro de suavização

Variando o valor de h obtemos diferentes formas de fh(x). Nos extremos,digamos, quando h → 0, temos uma representação muito ruidosa dosdados. Na situação oposta, quando h → ∞, temos uma representaçãomuito suave dos dados

A idéia do histograma serve como base para um estimador de densidadesmais geral conhecido como estimador naive (Silverman 1986). Seja Xuma v.a. com densidade f . Então,

f(x) = limh→0

1

2hP (x − h < X < x + h)

Para h fixo, podemos estimar P (x − h < X < x + h) pela proporção deobservações da amostra pertencentes ao intervalo (x − h, x + h). Dessemodo, um estimador natural de f , escolhendo h pequeno, é

f(x) =1

2nh#{Xi ∈ (x − h, x + h)}




Para expressar este estimador de forma mais clara, seja a função peso w:

w(x) =

{12

se |x| < 10 caso contrário.

(1)

Então, é fácil ver que uma estimativa para f neste caso é dada por

f(x) =1

n

n∑

i=1

1

hw

(x − Xi

h

)(2)

A partir de (1) podemos notar que o estimador (2) é construídocolocando-se um retângulo de largura 2h e altura (2nh)−1 em cadaobservação e então somando para obter a estimativa f

Não é difícil notar que f não é uma função contínua e tem derivada nulaem todos os pontos exceto nos pontos de salto X ± h

O estimador de densidades baseado em uma função kernel é obtidosubstituindo a função peso w por uma função não-negativa k, denominadafunção kernel, satisfazendo a condição

∫∞−∞ K(x)dx = 1

Usualmente, mas não sempre, K será uma função densidade deprobabilidade simétrica (Por exemplo, a função densidade de probabilidadenormal)




Estimação de Densidades Univaridas pelo Método Kernel

No caso univariado o estimador kernel para uma amostra aleatóriaX1, . . . , Xn retirada de uma distribuição com densidade comum f , podeser definido como

f(x; h) =1

nh

n∑

i=1

K

(x − Xi

h

)=

1

n

n∑

i=1

Kh (x − Xi) , (3)

onde h é o parâmetro de suavização, positivo e não-aleatório, e K é afunção kernel, não-negativa, satisfazendo a condição

∫ +∞−∞ K(x)dx = 1

A relação entre K e Kh é dada por Kh(t) = h−1K(h−1t)

Em cada ponto, uma função kernel dimensionada Kh com massa deprobabilidade n−1 é colocada. Estas são então somadas para fornecer acurva composta

A escolha da função kernel não é crucial para a performance do método, eé mais razoável escolher um kernel que auxilie na eficiência computacional(Silverman 1986; Epanechnikov 1969)




Tabela: Funções kernel comumente utilizadas com dados univariados

Função kernel Forma analítica, K(x)

Retangular 12

para |x| < 1, 0 caso contrário

Triangular 1 − |x| para |x| < 1, 0 caso contrário

Biweight 1516

(1 − x2)2 para |x| < 1, 0 caso contrário

Normal 1√2π

exp(−x2

2

)

Epanechnikov 34

(1 − x2/5

)/√

5 para |x| <√

5, 0 caso contrário




Figura: Estimativa da densidade univariada pelo método kernel. Linha sólida:densidade estimada; Linhas tracejadas: funções kernel individuais. A amostra écomposta pelos valores X1 = −1.0, X2 = −0.8, X3 = −0.6, X4 = 0.5, X5 = 1.2.Função kernel: gaussiana

−2 −1 0 1 2

0.0

0.1

0.2

0.3

0.4

0.5

0.6

x

dens

idad

e es

timad

a




Estimação de Densidades Multivaridas pelo Método Kernel

A extensão para dados multivariados é direta, com o estimador dedensidades p-dimensional, para uma amostra aleatória X

˜ 1, X˜ 2, . . . , X˜ n

retirada de uma densidade comum f , definido por

f(x˜) =

1

nhp

n∑

i=1

K

(1

h(x˜− X˜ i)

), (4)

onde x˜

= (x1, x2, . . . , xp)′ e X˜ i = (Xi1, Xi2, . . . , Xip)′, i = 1, 2, . . . , n

A função kernel multivariada K(x˜) é agora uma função definida no espaço

p-dimensional, satisfazendo∫Rp K(x

˜)dx˜

= 1

Usualmente K será uma função densidade de probabilidade unimodalradialmente simétrica




Exemplos de funções kernel multivariadas são a distribuição normal padrãomultivariada

K(x˜) = (2π)−p/2 exp

(−1

2x˜′x˜

),

e a função kernel Bartlett-Epanechnikov

K(x˜) =

{(1−x˜

′x˜)(p+2)

2cppara |x˜| < 1

0 caso contrário,

onde

cp =πp/2

Γ((p/2) + 1)

é o volume de uma esfera unitária p-dimensional




O uso de um único parâmetro de suavização em (4) implica que a funçãokernel colocada em cada ponto é dimensionada igualmente em todas asdireções e isso pode ser inadequado em muitas situações

Uma forma da estimativa da função de densidade de probabilidadecomumente utilizada é a soma do produto de funções kernel (sem,contudo, a implicação de independência entre as variáveis)

f(x˜) =

1

n

1

h1 · · ·hp

n∑

i=1

p∏

j=1

Kj

(xj − Xij

hj

), (5)

onde existem diferentes parâmetros de suavização associados com cadavariável. Pode-se assumir algum kernel univariado para os Kj ,j = 1, . . . , p. Usualmente, a mesma forma é assumida para todos os Kj .




Uma forma geral para o estimador de densidades multivariado, para umaamostra aleatória X

˜ 1, X˜ 2, . . . , X˜ n, retirada de uma densidade comum f , édada por

f(x˜) = f(x

˜;H) =

1

n

n∑

i=1

KH(x˜− Xi

˜), (6)

onde KH = |H|−1/2K(H−1/2x˜) é a função kernel dimensionada e H é

uma matrix não-aleatória, simétrica, positiva-definida, denominada matrizsuavização

A idéia básica do caso univariado, de colocar uma função kernel commassa de probabilidade n−1, é também válida no caso multivariado




Figura: Estimativa da densidade bivariada pelo método kernel. Linha sólida: curvas denível da densidade estimada; Linhas tracejadas: curvas de nível das funções kernelindividuais; Amostra: X

˜ 1 = (7, 3), X˜ 2 = (2, 4), X

˜ 3 = (4, 4), X˜ 4 = (5, 2),

X˜ 5 = (5.5, 6.5); Função kernel: gaussiana; Matriz de suavização: H =

[1 0.7

0.7 1

]

x

y

0 2 4 6 8 10

02

46

810

x

y

0 2 4 6 8 10

02

46

810




Escolha do Parâmetro de Suavização

O critério de erro mais amplamente utilizado nesta área de pesquisa é o(Erro Quadrático Integrado Médio) (EQIM) (Rosenblatt 1956), definidocomo

EQIM(f(·;H)) = E

{∫

Rp

[f(x˜;H) − f(x

˜)]2dx

˜

}(7)

Nosso objetivo é encontrar H tal que o EQIM seja minimizado, ou seja,

HEQIM = arg minH∈H

EQIM(f(·;H)), (8)

onde H é o espaço das matrizes simétricas, positivas-definidas dedimensão (p × p)

Contudo, o EQIM apresenta forma fechada apenas se f é uma mistura dedistribuições normais e K é a função kernel normal, e dessa forma,encontrar HEQIM é, em geral, extremamente difícil(Wand and Jones 1995)




O Erro Quadrático Integrado Médio Assintótico (EQIMA) é umaaproximação assintótica do EQIM. Uma expressão para o EQIMA,derivada por (Wand and Jones 1995), é

EQIMA(f(·;H)) =1

nR(K)|H|−1/2 +

1

4µ2(K)2

∫

Rp

tr2(HD2f(x˜))dx˜,

(9)onde R(v) =

∫Rp v(x

˜)2 para alguma função integrável quadrada v;

µ2(K)Ip =∫Rp x˜x˜′K(x

˜), com µ2(K) < ∞ e Ip é a matriz identidade de

dimensão (p × p); e D2f(x˜) é a matrix Hessiana de f

Devemos então encontrar um estimador do EQIM(A), EQIM(A), a partirdos dados disponíveis e a partir desse estimador encontrar um parâmetrode suavização H tal que

H = arg minH∈H

EQIM(A) (10)




A expressão (10) é chamada de seletor do parâmetro de suavização

Existem diversas metodologias que podem ser utilizadas para selecionar oparâmetro de suavização através de (10), dentre as quais, um métodoconhecido como plug-in e o método de mínimos quadrados por validaçãocruzada (LSCV, sigla em inglês)

Um estudo detalhado sobre métodos de seleção do parâmatro desuavização pode ser encontrado em (Duong 2004)




Análise Discriminante Não-paramétrica

De acordo com a regra de Bayes, nós alocamos uma observação para aclasse com maior probabilidade a posteriori:

x˜

é alocado para a classe Πj se Πj = arg maxj∈{1,...,J}

πjfj(x˜).

As probabilidades a priori πj , quando desconhecidas, podem ser estimadasusando πj = nj/n, j = 1, . . . , J , com

∑Jj=1 nj = n

Na metodologia paramétrica são feitas suposições sobre as densidades fj .

Usualmente, supõe-se que os dados seguem distribuição normal,entretanto, esta suposição pode ser muito restritiva ou até mesmoinadequada

Na análise discriminante não-paramétrica nós relaxamos essa suposiçãopara, dessa forma, poder lidar com casos mais complexos




−2 0 2 4

−2

02

4

x

y




−5 0 5 10

−5

05

x

y




A abordagem kernel para análise discriminante é estimar a densidade fj decada classe Πj e alocar uma observação de acordo com a regra:

x˜

é alocado para a classe Πj se Πj = arg maxj∈{1,...,J}

πj fj(x˜),

onde fj(x˜) é a estimativa da densidade pelo método kernel correspondente

a j-ésima classe




−2 0 2 4

−2

02

4

x

y

25

25

50

50 75

75

25

25

50

50

75

75




−5 0 5 10

−5

05

x

y




Suporte Computacional: R







Gratuito (disponível em http://www.R-project.org)

Código aberto

ColaborativoCentenas de pacotes implementados

AER: Applied Econometrics with RAMORE: A MORE flexibly neural networks packageAdMit: Adaptive Mixture of Student-t distributionsarules: Mining Association Rules and Frequent Itemsetsanapuce: Tools for microarray data analysisbetareg: Beta Regressionboot: Bootstrap R FunctionsBayesTree: Bayesian Methods for Tree Based Modelsclass: Functions for ClassificationclusterGeneration: Random cluster generation (with specified degree ofseparation)experiment: R package for designing and analyzing randomizedexperimentsFactoMiner: Factor Analysis and Data Mining with Rforeign: Read Data Stored by Minitab, S, SAS, SPSS, Stata, Systat,dBase, ...geoR: Analysis of geostatistical data




HiddenMarkov: Hidden Markov Models

intervals: Tools for working with points and intervals

JGR: Java Gui for R

kernlab: Kernel-based Machine Learning Lab

ks: Kernel density estimate for multivariate data

lodplot: Plot a genome scan

mcmc: Markov Chain Monte Carlo

nnet: Feed-forward Neural Networks and Multinomial Log-Linear Models

outliers: Tests for outliers

polspline: Polynomial spline routines

qcc: Quality Control Charts

ROCR: Visualizing the performance of scoring classifiers

survival: Survival analysis, including penalised likelihood

tree: Classification and regression trees

urca: Unit root and cointegration tests for time series data

VaR: Value at Risk estimation

e a lista cresce a cada dia...




ks: Kernel density estimate for multivariate data

> ## bivariate example

> data(unicef)

> H.scv <- Hscv(x=unicef)

> fhat <- kde(x=unicef, H=H.scv)

> plot(fhat, drawpoints=TRUE, drawlabels=FALSE, col=3, lwd=2)

50 100 150 200 250 300

4045

5055

6065

70

Under−5

Ave

life

exp




> plot(fhat, display="persp", border=NA, col="grey96",

+ shade=0.75)

Under−5

Ave life exp

Density function




> plot(fhat, display="image", col=rev(heat.colors(100)))

−100 0 100 200 300 400

3040

5060

7080

Under−5

Ave

life

exp




Uma função particularmente útil do pacote ks é a rmvnorm.mixt com a qualpodemos gerar dados oriundos de misturas de distribuições gaussianasmultivariadas.> mus <- rbind(c(-3/2,0), c(3/2,0))

> Sigmas <- rbind(diag(c(1/16, 1)), rbind(c(1/16, 1/18), c(1/18,

1/16)))

> props <- c(2/3, 1/3)

> x <- rmvnorm.mixt(1000, mus, Sigmas, props)

> plot(x, xlab = "x", ylab = "y")




−2 −1 0 1 2

−4

−3

−2

−1

01

2

x

y




> mus <- rbind(c(-2,3), c(2,3), c(0,2), c(0,1/2))

> Sigmas <- rbind(diag(c(1/10,1/10)), diag(c(1/10,1/10)),

diag(c(1/32,1/16)), diag(c(1,1/64)))

> props <- c(1/4, 1/4, 1/4, 1/4)

> x <- rmvnorm.mixt(1000, mus, Sigmas, props)

> plot(x, xlab = "x", ylab = "y")




−3 −2 −1 0 1 2 3

01

23

x

y




Experimentos Numéricos

Foram gerados dados simulando problemas com duas classes em seiscenários distintos

Amostras de treinamento de tamanhos 50, 100, 500 e 1000, gerados apartir de distribuições normais ou de misturas de distribuições normais

Amostras de teste independentes, fixas, de tamanho 1000

1000 réplicas de Monte Carlo

Foram comparados os métodos discriminante linear, quadrático e kernel

Métrica de comparação: taxa de erro no conjunto de teste

Todos os resultados foram obtidos através da linguagem R




Voltando...

Cenários Distribuição

A Π1 : f1 ∼ N

([1

−1

];

[49

1445

1445

49

]); Π2 : f2 ∼ N

([−1

1

];

[49

0

0 49

])

B Π1 : f1 ∼ N

([−1

1

];

[23

15

15

13

]); Π2 : f2 ∼ N

([1

−12

];

[23

15

15

13

])

C Π1 : f1 ∼ N

([−112

];

[23

15

15

49

]); Π2 : f2 ∼ N

([1

−12

];

[23

15

15

49

])

DΠ1 : f1 ∼

12

N

([−

32

−32

];

[45

−12

−12

45

])+ 1

2N

([1212

];

[45

−12

−12

45

]);

Π2 : f2 ∼12

N

([3232

];

[45

−12

−12

45

])+ 1

2N

([−

12

−12

];

[45

−12

−12

45

])

EΠ1 : f1 ∼

12

N

([−

32

0

];

[112

14

14

1

])+ 1

2N

([320

];

[112

14

14

1

]);

Π2 : f2 ∼ N

([0

0

];

[49

15

15

49

])

FΠ1 : f1 ∼

12

N

([−

32

0

];

[210

14

14

310

])+ 1

2N

([320

];

[310

14

14

310

]);

Π2 : f2 ∼ N

([0

0

];

[45

25

25

1

])




A

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

B

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

C

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

D

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

E

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

F

X1

X2

−3 −2 −1 0 1 2 3

−3

−2

−1

01

23




−3 −2 −1 0 1 2 3

−3

−2

−1

01

23

4

X1

X2

Π1Π2

−6 −4 −2 0 2 4

−4

−2

02

46

X1

X2

Π1Π2




Tabela: Resultados para os cenários A, B e C

Cenário Método Estatísticas n = 50 n = 100 n = 500 n = 1000

A

LDAMédia 0,6233 0,5664 0,5033 0,5039

D. P. 0,3017 0,2440 0,1081 0,0717

QDAMédia 0,3746 0,3381 0,3117 0,3025

D. P. 0,1790 0,0926 0,0451 0,0206

KDAMédia 1,8073 0,6713 0,3413 0,3244

D. P. 3,6037 1,0796 0,0962 0,0655

B

LDAMédia 1,0867 0,9642 0,9050 0,9088

D. P. 0,2814 0,1514 0,0736 0,0706

QDAMédia 1,3149 1,0702 0,9087 0,9052

D. P. 0,4757 0,2514 0,0746 0,0687

KDAMédia 1,5593 1,2551 0,9925 0,9410

D. P. 0,6088 0,3606 0,1576 0,1185

C

LDAMédia 4,8170 4,6447 4,5223 4,5146

D. P. 0,4767 0,3727 0,2348 0,1913

QDAMédia 4,9875 4,7300 4,5377 4,5195

D. P. 0,6089 0,4236 0,2559 0,2101

KDAMédia 2,5351 2,2039 1,8838 1,8341

D. P. 1,0339 0,7230 0,3798 0,2880




Tabela: Resultados para os cenários D, E e F

Cenário Método Estatísticas n = 50 n = 100 n = 500 n = 1000

D

LDAMédia 43.3035 44.5808 46.3690 46.6073

D. P. 2.6130 2.0645 0.5990 0.3108

QDAMédia 43.1092 44.0681 46.0483 46.4462

D. P. 2.6529 1.9639 0.7543 0.4151

KDAMédia 24.7502 19.9339 16.3670 16.1180

D. P. 6.4371 3.3978 0.6678 0.5049

E

LDAMédia 47.8123 48.3086 49.2980 49.5784

D. P. 3.2952 2.4899 1.4873 1.2574

QDAMédia 12.3440 10.0298 7.9060 7.7968

D. P. 4.0999 2.4674 0.6965 0.4928

KDAMédia 10.4084 7.5744 5.9294 5.9421

D. P. 3.8523 2.0576 0.5980 0.4642

F

LDAMédia 48.8693 49.1778 49.6406 49.6833

D. P. 2.4543 1.9858 1.1794 1.0829

QDAMédia 23.3331 21.6190 20.0978 19.8070

D. P. 4.6566 3.7946 2.1063 1.6127

KDAMédia 17.8189 15.4789 13.5814 13.1674

D. P. 2.8313 1.5302 0.5391 0.3796




Conclusões

Desempenho superior do método kernel em situações complexas (cenáriosD, E e F)

Desempenho similar aos métodos linear e quadrático em situações debaixa complexidade (cenários A, B e C)

Teoria bastante desenvolvida e consolidada

Implementações em linguagem R




Referências Bibliográficas

Duong, T. (2004).Bandwidth selectors for multivariate kernel density estimation.Ph. D. thesis, University of Western Australia, School of Mathematics and Statistics.

Epanechnikov, V. A. (1969).Non-parametric estimation of a multivariate probability density.Theory of Probability and its Applications 14, 153–158.

Rosenblatt, M. (1956).Remarks on some nonparametrics estimates of a density function.The Annals of Mathematical Statistics. 27, 832–837.

Scott, D. W. (1992).Multivariate Density Estimation: Theory, Practice, and Visualization.New York: John Wiley & Sons.

Silverman, B. W. (1986).Density Estimation for Statistics and Data Analysis.London: Chapman & Hall.

Simonoff, J. S. (1996).Smoothing Methods in Statistics.New York: Springer-Verlag.

Wand, M. P. and M. C. Jones (1995).Kernel Smoothing.London: Chapman & Hall.



Documents

M todo Kernel: Estima o de Densidades e Classifica o de ...fatc/AM/kernel.pdf · Kernel (1) NÃO é o núcleo celular NÃO é o núcleo de um sistema operacional NÃO é o núcleo/espaço