Estimação dos parâmetros da mistura de duas componentes GEV via Algoritmo EM

Estimação dos parâmetros da mistura de duas componentes GEV via Algoritmo EM^*

C.E.G. Otiniano^**; E.C.M. Teixeira

Departamento de Estatística, Universidade de Brasília, 70910-900 Brasília, DF, Brasil. E-mail: [e-mail]

RESUMO

Modelos probabilísticos de misturas finitas de densidades cujas componentes são as densidades da distribuição de Valor Extremal Generalizado tem uma ampla aplicabilidade em diversas áreas como finanças e hidrologia. Neste trabalho, os estimadores dos parâmetros da mistura de duas componentes GEV são obtidos via o algoritmo EM. Apresentamos também ilustrações numéricas do comportamento dos estimadores obtidos através de simulação, bem como uma aplicação a dados reias.

Palavras-chave: GEV 1, mistura finita 2, algoritmo EM 3.

ABSTRACT

Probabilistic models of finite mixtures with components Generalized Extremal Value (GEV) distributions are widely applied in diverse areas such as finance and hydrology. In this work, the estimators of the parameters of the GEV mixture of two components are obtained via the EM algorithm. We also present numerical examples of the behavior of the estimates obtained through of simulation, well as an application to real data.

Keywords: GEV, finite mixture, EM algorithm.

1 INTRODUÇÃO

Dadas as funções de densidade de probabilidade (f.d.p.) d-dimensionais f₁, ..., f_k de uma família de densidades

= {f (x; θ), x ∈ x R^d, θ ∈ Θ},

onde Θ é o espaço paramétrico dessa família. Uma mistura finita de f₁, ..., f_k ∈ é uma f.d.p.

onde os números p₁, ..., p_k são chamados de pesos da mistura e f₁, ..., f_k chamadas decomponentes da mistura.

Devido a sua definição, as misturas finitas de distribuições ou de densidades de probabilidadesão frequentemente utilizadas para modelar dados de populações heterogêneas. Tais misturas permitem construir modelos probabilísticos de uma ampla variedade de fenômenos, em por exemplo, engenharia, economia e hidrologia, entre outros. Teoria e aplicações de misturas podem ser encontradas em Titterington et al. (1985), McLachlan & Basford (1998), e McLachlan & Peel (2000).

Neste trabalho abordamos a estimação do modelo (1.1) quando d = 1 e k = 2, onde as componentes da mistura são densidades da distribuição de Valor Extremal Generalizado univariadas.

A distribuição de Valor Extremal Generalizado do inglês generalized extreme value (GEV) distribution, G_γ, é a distribuição limite de máximos (ou mínimos) normalizados de sequências de variáveis aleatórias (v.a.'s) independentes e identicamente distribuídas (i.i.d.). Isto é, seja {X_n} tal sequência. Se existirem sequências de números reais {c_n}, c_n > 0 e {d_n} tais que

então a distribuição GEV introduzida por Jenkinson (1955) é definida por

onde γ, µ ∈ R e σ > 0 são parâmetros de forma (índice caudal), locação e escala, respectivamente. Aqui 1 + (x - µ) > 0.

A correspondente f.d.p. de G_γ, é dada por

Estas distribuições são utilizadas em hidrologia para analisar a frequência de fluídos e em finanças para calcular o valor em risco (VaR) de retornos, perdas ou ganhos extremos.

Escalante-Sandoval (2007) utilizou misturas de duas componentes GEV para avaliar dados sobre a frequência de fluídos de rios de uma região do Noroeste do México. Neste caso, os parâmetros foram estimados por máxima verossimilhança utilizando vários métodos computacionais intensivos e as estimativas obtidas por esse procedimento, para vários casos, não convergiram. Outro trabalho que envolve mistura de distribuições com alguma componente GEV é o de Kollu et al. (2012), no qual, com o objetivo de descrever características da velocidade de ventos, os autores apresentam modelos de mistura de duas densidades cujas componentes são Weibull e GEV, Weibull e lognormal, e GEV e lognormal. Neste caso, a estimação também foi feita por máxima verossimilhança. Em geral, o método de máxima verossimilhança, o algoritmo EM, e os métodos bayesianos são os métodos alternativos de estimação dos parâmetros de misturas finitas de distribuições.

Com o objetivo de apresentar um método de estimação alternativo ao de Escalante-Sandoval (2007), dos parâmetros de misturas de duas componentes GEV, neste trabalho descrevemos as expressões a serem calculadas utilizando o algoritmo EM e mostramos o resultado das estimativas para alguns experimentos. Uma aplicação com dados reais também é adicionada.

Este trabalho esta organizado em cinco seções. Na Seção 2, apresentamos os principais resultados deste trabalho, expressões que devem ser calculadas para estimar os parâmetros da mistura de duas componentes GEV via algoritmo EM. Na seção 3, mostramos os resultados das estimações obtidas para algumas simulações. Já na seção 4, uma aplicação da estimação de um modelo de mistura para dados reais é apresentada. Por fim, nas Seções 5 e 6 fazemos os agradecimentos e citamos as referência bibliográficas, respectivamente.

2 ESTIMAÇÃO VIA ALGORITMO EM

O algoritmo EM do inglês "Estimation Maximization" é um algoritmo clássico da estatística usado para determinar estimativas de parâmetros em modelos de mistura, em modelos com dados faltantes, e em modelos em que a estimação dos parâmetros por máxima verossimilhança apresenta problemas. O livro McLachlan & Peel (2000) é uma boa referência para esse assunto. Nesta seção utilizamos o algoritmo EM para obter as estimativas de misturas de duas componentes GEV.

Seja X uma v.a. com f.d.p. h(x, Θ) mistura de duas componentes, definida em (1.1), e dada por

onde Θ = (p₁, θ₁, θ₂) com θ_ℓ = (γ_ℓ, σ_ℓ, µ_ℓ), ℓ = 1, 2 e γ_ℓ(.; θ_ℓ) é uma componente de uma das famílias

Nesta seção, usamos o Algoritmo EM para obter estimativas de Θ, baseado nos valores x₁, x₂, ..., x_N de uma amostra aleatoria de tamanho N da v.a. X.

Uma breve descrição do algoritmo EM, é a seguinte. Se assumirmos que X é observado e gerado por alguma distribuição paramétrica f(x; Θ), chamamos X ({x₁, ..., x_N}) de dados incompletos. Estes dados são completados com Y({y₁, ..., y_N}), então Z = (X, Y) são os dados completos cuja densidade conjunta é

f(z; Θ) = f(x, y; Θ) = f(y / x, Θ) f(x; Θ).

Com esta nova densidade, defini-se a função de verossimilhança dos dados completos

L(Θ / z) = L(Θ / x, y) = f(X, Y; Θ).

O algoritmo EM alterna o Passo E (da esperança), onde obtem-se

Q(Θ, Θ^(k)) = E[1n(f(X, Y / Θ)) / X, Θ^(k)],

com o Passo M (da maximização), onde se calcula Θ^{(k + 1)} ao maximizar Q(Θ, Θ^(k)).

No caso de mistura, f(x; Θ) = h(x; Θ) dada em (2.1) e Q é dada pela expressão

Em síntese, na (k + 1)-ésima iteração do passo-E, a atualização da estimativa de p₁, é dada por

onde,

e na (k + 1)-ésima iteração do passo-M, a atualização das estimativas de , ℓ = 1,2 são obtidas resolvendo as equações

Na estimação de Θ do modelo (2.1), como cada componente da mistura g_ℓ(.; θ_ℓ) deve ser de uma das famílias, G_γ+, G_γ-, G₀, há seis possíveis casos para aplicar. Os três primeiros casos descritos a seguir correspondem a mistura de componentes da mesma família e os outros três casos seguintes correspondem a mistura de componentes das diferentes famílias.

Caso 1. (Mistura de componentes em ₀) As componentes g₁(.; θ₁) e g₂(.; θ₂) no modelo (2.1) são da família ₀, onde θ_ℓ = (σ_ℓ, µ_ℓ), ℓ = 1,2, e Θ = (p₁, θ₁, θ₂).

No passo-E, as atualizações de são obtidas conforme (2.2). As atualizações de , ℓ = 1,2 são obtidas ao resolver, utilizando o método de Newton-Raphson, a seguinte equação:

Para as atualizações de , ℓ = 1,2 usamos a fórmula fechada

Caso 2. (Mistura de componentes em ₊). As componentes g₁(.; θ₁) e g₂(.; θ₂) no modelo (2.1) são da família G₊, onde θ_ℓ = (γ_ℓ, σ_ℓ, µ_ℓ), ℓ = 1,2, e Θ = (p₁, θ₁, θ₂).

Quando µ₁ - = µ₂ - , no passo-E, as atualizações de são obtidas conforme (2.2) e (2.3). Porém, quando µ₁ - < µ₂ - , as atualizações de são obtidas de (2.2), onde (2.3) é substituido por

e

No passo-M, as atualizações de , , e , ℓ = 1,2 são obtidas ao resolver as equações

e

respectivamente.

Caso 3. (Mistura de componentes em _–). As componentes g₁(.; θ₁) e g₂(.; θ₂) no modelo (2.1) são da família _–, onde θ_ℓ = (γ_ℓ, σ_ℓ, µ_ℓ), ℓ = 1,2, e Θ = (p₁, θ₁, θ₂).

Quando µ₁- = µ₂-, no passo-E, as atualizações de são obtidas conforme (2.2) e (2.3). Porém, quando µ₁-< µ₂-, em vez de (2.3) utilizamos

e

No passo-M, as atualizações de , , e , ℓ = 1,2 são obtidas ao resolver as equações (2.9), (2.10), e (2.11), respectivamente.

Caso 4. (Mistura de componentes em ₊∪ ₀ ). As componentes no modelo (2.1) são g₁(x; θ₁) ∈ ₊ e g₂(x; θ₂) ∈ ₀, onde θ₁ = (γ₁, σ₁, µ₁), e θ₁ = (σ₂, µ₂) ℓ = 1,2. Então Θ = (p₁, γ₁, σ₁, µ₁, σ₂, µ₂).

No passo-E, as atualizações de são obtidas conforme (2.2), onde

e

No passo-M, as atualizações de , , e são obtidas ao resolver as equações (2.9), (2.10), e (2.11), respectivamente. Para as atualizações de e utilizamos (2.5) e (2.6), respectivamente.

Caso 5. (Mistura de componentes em _– ∪ ₀). As componentes no modelo (2.1) são g₁(x; θ₁) ∈ _– e g₂(x; θ₂) ∈ ₀, onde θ₁ = (γ₁, σ₁, µ₁), e θ₁ = (σ₂, µ₂) ℓ = 1,2. Então Θ = (p₁, γ₁, σ₁, µ₁, σ₂, µ₂).

No passo-E, as atualizações de são obtidas conforme (2.2), onde

e

No passo-M, as atualizações de , , , e são obtidas como no Caso 4.

Caso 6. (Mistura de componentes em _–∪ ₊). Neste caso, g₁(x; θ₁) ∈ ₊ e g₂(x; θ₂) ∈ _– onde θ_ℓ = (γ_ℓ, σ_ℓ, µ_ℓ), ℓ = 1,2 and Θ = (p₁, γ₁, σ₁, µ₁, γ₂, σ₂, µ₂).

Consideramos os dois possíveis sub casos:

Quando µ₁-< µ₂-

e

E quando µ₁-> µ₂-

and

Neste caso as atualizações de , são obtidas conforme , , e , ℓ = 1,2 são obtidas ao resolver as equações (2.2), (2.9), (2.10), e (2.11), respectivamente.

3 SIMULAÇÃO

Nesta seção, testamos via simulação as estimativas de Θ do modelo (2.1) para a mistura de componetes de diferentes famílias, o que corresponde trabalhar com os casos (4)-(6), pois os casos (1)-(3) tratam de componentes da mesma família que são mais facilmente calculados. Para isto, seguimos o seguinte procedimento:

1. Geramos amostras aleatórias de tamanho n = 100 para cada escolha do vetor Θ.

2. A amostra aleatória da variável aleatória X cuja densidade é a mistura (2.1) é gerada como segue:

(a) Gerar duas variáveis uniformes u1 e u2.

(b) Se u₁ < p₁, então usamos u₂ para gerar um valor x da v.a. X de (2.1), onde x = (u₂) e G₁ é a distribuição acumulada g₁.

(c) Se u₁> p₁, então usamos u₂ para gerar um valor x da v.a. X, onde x = (u₂) e G₂ é a distribuição acumulada g₂.

3. Calculamos iterativamente os estimadores de Θ utilizando as expressões de g(ℓ / x_i, θ^(k)), ℓ = 1,2 mostradas nos casos (4)-(6).

4. A regra de parada que utilizamos no Algoritmo EM é log(L[Θ^{(k + 1)}]) - log(L[Θ^(k)]) < n 10^-5.

5. Obtemos estimativas de Θ utilizando 100 amostras de tamanho 100, então calculamos a média e o erro quadrático médio de Θ por Monte Carlo. Tais resultados são mostrados nas Tabelas 2 e 3.

Definimos vários experimentos, valores de Θ, para testar o comportamento dos estimadores do algoritmo EM. Esses experimentos mostrados na Tabela 1 foram escolhidos de tal maneira a contemplar os casos (4)-(6).

Os experimentos Θ_4,i, Θ_5,i, Θ_6,ia, e Θ_6,ib, i = 1, 2, 3, correspondem aos casos (4), (5), e (6), respectivamente.

Nas Figuras 1-4, comparamos a densidades de uma amostra h(x; Θ) com a densidade h(x; ), onde os valores de são as médias de , mostradas na Tabela 2. Por essas figuras podemos concluir que o algoritmo EM teve um bom desempenho, pois a densidade h(x; ) teve um bom ajuste, em quase todos os experimentos exceto no caso Θ_6,2b. Os resultados do erro quadrático médio das estimativas são mostrados na Tabela 3, os quais confirmam que os resultados da Tabela 1 são bons estimadores dos parâmetros em análise. Novamente apenas o caso Θ_6,2b não foi muito bom.

4 APLICAÇÃO

Para demonstrar a aplicabilidade da mistura (2.1), ajustamos o logaritmo do consumo per capita de petróleo em 135 países no ano de 2001. Os dados foram obtidos no endereço http://data.worldbank.org/indicator/EG.USE.PCAP.KG.OE.

Para esses dados, os valores dos parâmetros estimados do modelo (2.1), pelo algoritmo EM, foram p = 0.5, σ₁ = 0.5155666, σ₂ = 0.5985117, µ₁ = 6.414, µ₂ = 8.235. O que indica a mistura de duas componentes da família ₀. Na Figura 5 mostra-se a densidade ajustada aos dados.

O QQ-Plot da Figura 6, mostra que o ajustamento da mistura de duas distribuições com densidade Gumbel é adequado.

5 CONCLUSÕES

Neste artigo, os estimadores dos parâmetros do modelo (2.1) são dados em seis possíveis casos, usando o Algoritmo EM. O método de estimação é testado com doze experimentos. Apenas para um desses experimentos o erro quadrático médio não foi suficientemente pequeno. Isso pode ter ocorrido pelo fato dos parâmetros de locação estarem muito distantes. Finalmente, utilizamos dados reais para mostrar uma aplicação do algoritmo EM. Não fizemos comparação de nossas estimativas com as de por exemplo Atienza-Sandoval (2007), pois os algoritmos por eles utilizados demandam um custo computacional intensivo. O tempo computacional para o cálculo de nossas estimativas foi inferior a um minuto.

AGRADECIMENTOS

Agradecemos aos revisores deste artigo por seus valiosos comentários para melhorar a primeira versão. Agradecemos também as agências de fomento à pesquisa CAPES/PROCAD e DPP/UnB pelo financiamento parcial para este trabalho.

Recebido em 9 julho, 2013

Aceito em 26 março, 2014

• [1] C. Escalante-Sandoval. A Mixed distribution with EV1 and GEV components for analyzing heterogeneous samples. Ingeniería Investigación y Tecnología 8, 3 (2007), 123-133.
• [2] R. Kollu, S.R. Rayapudi, S. Narasimham & K.M. Pakkurthi. Mixture probability distribution functions to model wind speed distributions. International Journal of Energy and Environmental Engineering, (2012), 1-10.
• [3] G.J. McLachlan & K.E. Basford. "Mixture models". Marcel Dekker, NewYork, (1988).
• [4] G.J. McLachlan & D. Peel. "Finite mixture models". Wiley, New York, (2000).
• [5] D.M. Titterington, A.M.F. Smith & U.E. Makov. "Statistical analysis of finite mixture distributions". Wiley, NewYork, (1985).

*

Pesquisa parcialmente financiada por CAPES/PROCAD.

**

Autor correspondente: Cira Etheowalda Guevara Otiniano

Datas de Publicação

• Publicação nesta coleção
  10 Jun 2014
• Data do Fascículo
  Abr 2014

Histórico

• Recebido
  09 Jul 2013
• Aceito
  26 Mar 2014