Capítulo 3

Histogramas: teoria e prática


A estrutura do histograma clássico é útil para transmitir a essência da teoria e prática não paramétricas. Da teoria do erro quadrático aos algoritmos de validação cruzada de última geração, esses tópicos devem ser estudados cuidadosamente. Discutir essas questões no contexto do estimador de densidade não paramétrico mais simples e conhecido é de interesse prático e será valioso no desenvolvimento da teoria de estimadores mais complexos.

O histograma é geralmente apresentado em uma escala não-densitária: seja como contagens de classes ou como um diagrama de caule e folhas (Tukey 1977). Alguns autores estabeleceram uma distinção entre os propósitos de um histograma como estimador de densidade e como dispositivo de apresentação de dados (Emerson and Hoaglin 1983).

Tal distinção parece artificial, pois não há razão para desconsiderar um histograma construído de maneira menos que ideal. O exame de histogramas subsuavizados e supersuavizados deve ser rotineiro. Um histograma transmite informações visuais tanto da frequência quanto das frequências relativas das observações; ou seja, a essência de uma função de densidade.



3.1 Regra de Sturges para seleção da largura de banda


O histograma de frequência clássico é formado pela construção de um conjunto completo de intervalos não sobrepostos, chamados classes, e pela contagem do número de pontos em cada classe. Para que a contagem de itens nos recipientes seja comparável, todos os recipientes devem ter a mesma largura. Nesse caso, o histograma é completamente determinado por dois parâmetros: a largura do intervalo ou largura de banda \(h\) e a origem do intervalo \(t_0\), que é qualquer ponto final do intervalo escolhido convenientemente. Frequentemente, a origem do intervalo é escolhida como \(t_0 = 0\).

Embora a ideia de agrupar dados na forma de um histograma seja pelo menos tão antiga quanto o trabalho de Graunt em 1662, nenhuma diretriz sistemática para a construção de histogramas foi fornecida até a breve nota de Herbert Sturges em 1926. Seu trabalho utilizou um recurso que foi defendido de forma mais geral por Tukey (1977); ou seja, tomar a densidade normal como ponto de referência ao se pensar em dados. Sturges simplesmente observou que a distribuição binomial \(B(n, p = 0.5)\), poderia ser usada como um modelo de um histograma construído de forma otimizada com dados normais escalonados adequadamente (ver Figura 3.1).

Imaginemos construir um histograma de frequência com \(k\) classes, cada uma com largura 1 e centrada nos pontos \(i = 0, 1,\cdots, k-1\). Escolhemos a contagem de classes da \(i\)-ésima classe como a distribuição binomial \({k-1 \choose i}\). À medida que \(k\) aumenta, este histograma de frequência ideal assume a forma de uma distribuição normal com média \((k-1)/2\) e variância \((k-1)/4\).

O tamanho total da amostra é \[ n=\sum_{i=0}^{k-1} {k-1 \choose i} = (1+1)^{k-1} = 2^{k-1} \] pela expansão binomial. A regra de Sturges segue imediatamente: \[ \tag{3.1} \mbox{Regra de Sturges sobre o número de classes}:\quad k = 1+log_2(n)\cdot \]

É claro que qualquer múltiplo constante dos coeficientes binomiais também pareceria normal, mas observe que a condição de contorno \(n = 1\) resulta em \(k = 1\) na Equação (3.1).

Na prática, a regra de Sturges é aplicada dividindo-se o intervalo amostral dos dados no número prescrito de classes de largura igual. Tecnicamente, a regra de Sturges é uma regra de número de classes, e não de largura de classe. Uma abordagem muito mais simples é adotar a convenção de que todos os histogramas têm um número infinito de classes, das quais apenas um número finito não está vazio.

Além disso, histogramas adaptativos, que são considerados na Seção 3.2.8, não utilizam classes de largura igual. Portanto, o foco na largura da classe, em vez do número de classes, parece apropriado.

Figura 3.1: Função de densidade de probabilidade binomial com \(p = 0.5\) usada por Sturges para determinar o número de intervalos do histograma.

A regra de Sturges é amplamente recomendada em livros introdutórios de estatística e frequentemente utilizada como padrão em pacotes estatísticos. Se os dados não forem normais, mas assimétricos ou leptocurtóticos, podem ser necessários intervalos adicionais. Por exemplo, Doane (1976) propôs aumentar o número de intervalos em (3.1) para \(\log_2\left(1+\widehat{\gamma}\sqrt{n/6}\right)\), onde \(\widehat{\gamma}\) é uma estimativa do coeficiente de assimetria (skewness) padronizado (ver Exercício 3.1).


3.2 A teoria \(L_2\) de histogramas univariados


3.2.1 Erro quadrático médio pontual e consistência


Nesta seção, são apresentadas as propriedades do erro quadrático médio (MSE) de um histograma de densidade. A diferença entre um histograma de frequência e um histograma de densidade é que este último é normalizado para integrar igual a 1. Como mencionado anteriormente, o histograma é completamente determinado pela amostra \(\{x_1,\cdots,x_n\}\) de \(f(x)\) e por uma escolha de malha \(\{t_k,-\infty< k < \infty\}\). Seja \(B_k = [t_k,t_{k+1})\) o \(k\)-ésimo intervalo.

Suponha que \(t_{k+1}-t_k = h\) para todo \(k\); então, diz-se que o histograma tem largura de intervalo fixa \(h\). Um histograma de frequência é construído usando blocos de altura 1 e largura \(h\) empilhados nos intervalos apropriados. A integral de tal figura é claramente igual a \(nh\). Assim, um histograma de densidade usa blocos de construção de altura \(1/(nh)\), de modo que cada bloco tenha área igual a \(1/n\).

Seja \(\nu_k\) a contagem de amostras do \(k\)-ésimo bloco, ou seja, o número de pontos amostrais que caem no bloco \(B_k\) (veja Figura 3.2). Então, o histograma é definido da seguinte forma: \[ \tag{3.2} \widehat{f}(x)=\dfrac{\nu_k}{nh}=\dfrac{1}{nh}\sum_{i=1}^n \pmb{I}_{t_k,t_{k+1}}(x_i), \quad \mbox{para} \quad x\in B_k\cdot \]

A análise da variável aleatória do histograma \(\widehat{f}(x)\), é bastante simples, uma vez reconhecido que as contagens dos intervalos, \(\{\nu_k\}\), são variáveis aleatórias binomiais: \[ \nu_k \sim B(n,p_k), \quad \mbox{onde} \quad p_k = \int_{B_k} f(t)\mbox{d}t\cdot \]

Figura 3.2: Notação para a construção de um histograma com espaçamento uniforme.

Considere o MSE de \(\widehat{f}(x)\) para \(x\in B_k\). Agora \(\mbox{E}(\nu_k)=n \, p_k\) e \(\mbox{Var}(\nu_k)=n \, p_k(1-p_k)\). Então, \[ \tag{3.1} \mbox{Var}\big(\widehat{f}(x) \big) =\dfrac{\mbox{Var}(\nu_k)}{(n\, h)^2}=\dfrac{n \, p_k(1-p_k)}{(n\, h)^2} \] e \[ \tag{3.4} \mbox{Viés}\big(\widehat{f}(x) \big)=\mbox{E}\big(\widehat{f}(x) \big)-f(x)=\dfrac{1}{n\, h} \mbox{E}(\nu_k)-f(x)=\dfrac{p_k}{h}-f(x)\cdot \]

Para prosseguir com o mínimo de suposições, suponha que \(f(x)\) seja Lipschitz contínua no intervalo \(B_k\).


Definição 3.1: Diz-se que uma função \(f(x)\) é Lipschitz contínua em um intervalo \(B_k\) se existir uma constante positiva \(\gamma_k\) tal que \(|f(x)-f(y)|< \gamma_k |x-y|\) para todo \(x,y\in B_k\).


Então, pelo teorema do valor médio (MVT), \[ \tag{3.5} p_k = \int_{B_k} f(t)\mbox{d}t=h\, f(\xi_k), \qquad \mbox{para algúm} \qquad \xi_k\in B_k\cdot \] Conclui-se, portanto, que \[ \tag{3.6} \mbox{Var}\big( \widehat{f}(x)\big)\leq \dfrac{p_k}{n\, h^2}=\dfrac{f(\xi_k)}{n\, h} \] e \[ \left| \mbox{Viés}\big( \widehat{f}(x)\big) \right|=\left| f(\xi_k) -f(x)\right| \leq \gamma_k |\xi_k -x|\leq \gamma_k h; \] por isso, \(\mbox{Viés}^2\big( \widehat{f}(x)\big)\leq \gamma_k^2 \, h^2\) e \[ \tag{3.7} \mbox{MSE}\big( \widehat{f}(x)\big) \leq \dfrac{f(\xi_k)}{n\, h}+\gamma_k^2 \, h^2\cdot \]

Na literatura sobre suavização, a largura do intervalo \(h\) é referida como um parâmetro de suavização, pois controla o grau de suavidade no estimador para uma dada amostra de tamanho \(n\).

A equação (3.7) resume a relação inversa recorrente entre viés e variância, determinada pela escolha do parâmetro de suavização. A variância pode ser controlada tornando \(h\) grande, de modo que os intervalos sejam largos e de altura relativamente estável; no entanto, o viés é grande. Por outro lado, o viés pode ser reduzido tornando \(h\) pequeno, de modo que os intervalos sejam estreitos; no entanto, a variância é grande.

Observe que o viés pode ser eliminado escolhendo \(h = 0\), mas esse histograma muito grosseiro é exatamente a função de densidade de probabilidade empírica \(f_n(x)\), que tem variância infinita (vertical). O viés e a variância podem ser controlados simultaneamente escolhendo um valor intermediário para a largura do intervalo e permitindo que a largura do intervalo diminua lentamente à medida que o tamanho da amostra aumenta.


Definição 3.2: Diz-se que um estimador de densidade é consistente na média quadrática se \[ \lim_{n\to\infty} \mbox{MSE}\big( \widehat{f}(x)\big) = 0 \cdot \]


Um parâmetro de suavização ótimo \(h^*\) é definido como aquela escolha que minimiza o \(\mbox{MSE}\) (assintótico). Os resultados a seguir são consequências da equação (3.7).


Teorema 3.1: Suponha que \(x\in B_k\) seja um ponto fixo e que \(f\) seja Lipschitz contínua neste intervalo com \(\gamma_k\) constante. Então a estimativa do histograma \(\widehat{f}(x)\) é consistente em média quadrática se, quando \(n\to\infty\), então \(h\to 0\) e \(n h\to \infty\).


Demonstração Exercício.


A primeira condição garante que o viés desapareça assintoticamente, enquanto a segunda condição garante que a variância tenda a zero. Duda and Hart (1973) sugerem escolher \(h = n^{-1/2}\), por exemplo.


Corolário 3.2: O limite do \(\mbox{MSE}(x)\) (3.7) é minimizado quando \[ \tag{3.8} h^*(x)=\left( \dfrac{f(\xi_k)}{2\gamma_k^2\, n}\right)^{1/3}; \] o \(\mbox{MSE}^*(x)\) resultante é \(O(n^{-2/3})\).


Esses resultados merecem uma análise cuidadosa. A largura ideal do intervalo diminui a uma taxa proporcional a \(n^{-1/3}\). Essa taxa é muito mais rápida do que a regra de Sturges, que sugere a taxa \(\log^{-1}_2(n)\) (ver Tabela 3.1). A taxa ideal de diminuição do \(\mbox{MSE}\) não atinge o limite inferior de Cramer-Rao de \(O(n^{-1})\) para estimadores paramétricos.

O ruído inerente ao histograma varia diretamente com a raiz quadrada de sua altura, uma vez que \(\mbox{Var}\big(\widehat{f}(x)\big)\approx f(x)/(n\, h)\) da equação (3.6). Essa heterocedasticidade (variância desigual) na estimativa do histograma pode ser eliminada usando uma transformação de estabilização de variância. Cada contagem de intervalo é aproximadamente uma variável aleatória de Poisson. É bem conhecido que a raiz quadrada é a transformação de estabilização de variância para dados de Poisson.

Suponha que \(Y_n\) tenha momentos \((\mu_n,\sigma^2_n)\) com \(\mu_n > 0\) e \(\sigma_n\to 0\). Então \[ \mbox{Var}\big(g(Y_n)\big) \approx g'(\mu_n)^2 \sigma^2_n\cdot \] Escolhendo \(Y_n = \widehat{f}(x)\) de modo que \(\mu_n\approx f(x)\) e \(\sigma^2_n\approx f(x)/(n\, h)\), então \(g(y) =\sqrt{y}\) e \(g'(y) = 1/(2\sqrt{y})\) ou \(1/(2\sqrt{f(x)})\) em \(y = \mu_n\). Portanto, \[ \tag{3.9} \sqrt{\mbox{Var}\Big(\sqrt{\widehat{f}(x)} \Big) } \approx \dfrac{1}{2\sqrt{f(x)}}\sqrt{\dfrac{f(x)}{n\, h }} = \dfrac{1}{2\sqrt{n\, h}}, \] que é independente da incógnita \(f(x)\).

Assim, mostrar graficamente o histograma em escala de raiz quadrada permite uma comparação fácil do ruído no histograma em regiões de alta e baixa densidade. Tukey (1977) chamou a estimativa resultante de rootgrama. Obviamente, o rootgrama não retrata mais com precisão as frequências relativas das observações.

Em mais de uma dimensão, o rootgrama ainda estabiliza a variância. No entanto, como os contornos do histograma bivariado e do rootgrama bivariado são idênticos, as aplicações práticas são limitadas.

Uma consequência do Corolário 3.2 é que o uso de uma largura de banda fixa em toda a extensão dos dados geralmente não é ideal. Por exemplo, a largura do intervalo deve ser relativamente maior em regiões de maior densidade para reduzir a variância na equação (3.7). Agora, se a largura do intervalo \(B_k\) for suficientemente estreita, a constante de Lipschitz \(\gamma_k\) é essencialmente a magnitude da inclinação de \(f(x)\) nesse intervalo.

Portanto, a partir da equação (3.8), a largura do intervalo deve ser menor em regiões onde a densidade está mudando rapidamente, e vice-versa. Essas noções são confirmadas na Seção 3.2.8. No entanto, na prática, não existem algoritmos confiáveis para construir malhas de histograma adaptativas. Portanto, o estudo de histogramas de largura fixa continua sendo importante.


3.2.2 Erro global \(L_2\) do histograma


Resultados de consistência baseados em limites superiores não são úteis na prática, visto que esses limites podem estar bastante distantes da realidade. Aproximações mais úteis podem ser obtidas assumindo a existência de derivadas de \(f\). Esses resultados podem ser úteis mesmo em situações práticas onde \(f\) é desconhecido, empregando-se diversas técnicas chamadas algoritmos de validação cruzada (ver Seção 3.3).

O cálculo do erro quadrático médio integrado (\(\mbox{MISE}\)) é realizado agregando-se o MSE em cada intervalo e somando-se os valores de todos os intervalos. Considere-se primeiro a integral da variância (\(\mbox{IV}\)): \[ \tag{3.10} \mbox{IV} = \int_{-\infty}^\infty \mbox{Var}\big(\widehat{f}(x)\big)\mbox{d}x = \sum_{k=-\infty}^\infty \int_{B_k} \mbox{Var}\big(\widehat{f}(x)\big)\mbox{d}x\cdot \]

Da equação (3.3), a última integral sobre \(B_k\) é simplesmente \(p_k(1-p_k)/(nh)\). Agora \[ p_k = \displaystyle \int f(x)\mbox{d}x = 1\cdot \] Lembre-se que \[ \sum \phi(\epsilon_k)\times h = \int \phi(x)\mbox{d}x + o(1) \] pela aproximação integral Riemanniana padrão. Portanto, usando a aproximação (3.5) para \(p_k\), \[ \sum_k p_k^2 = \sum_k f(\epsilon_k)^2 h^2 = h \sum_k f(\epsilon_k)^2 h = h\left(\int f(x)^2\mbox{d}x+o(1) \right)\cdot \] Combinando, temos \[ \tag{3.11} \mbox{IV} = \dfrac{1}{nh}-\dfrac{R(f)}{n}+o\big(n^{-1}\big), \] onde a seguinte notação é adotada para a norma \(L_2\) ao quadrado de \(\phi\): \[ R(\phi)=\int \phi^2(x)\mbox{d}x\cdot \]

A norma $L_2 $ao quadrado é apenas uma das medidas possíveis da rugosidade \((R)\) da função \(\phi\). Alternativas incluem \(R(\phi ')\) e \(R(\phi '')\). \(R(\phi)\), neste contexto, refere-se mais à rugosidade estatística do que à rugosidade matemática da função \(\phi\).

Esta última geralmente se refere ao número de derivadas contínuas na função. A primeira se refere, de forma geral, ao número de oscilações na função. A rugosidade estatística, no entanto, leva em conta as primeiras derivadas da função. Em resumo, nenhuma das definições descreveria uma densidade normal como rugosa, mas a densidade lognormal é uma função muito rugosa do ponto de vista estatístico, mesmo sendo infinitamente diferenciável. Uma densidade polinomial de baixa ordem, como uma da família Beta, é estatisticamente suave, mesmo possuindo poucas derivadas.

Para calcular o viés, considere um intervalo típico \(B_0 = [0,h)\). A probabilidade do intervalo \(p_0\), pode ser aproximada por \[ \begin{array}{rcl} p_0 & = & \displaystyle \int_0^h f(t)\mbox{d}t = \int_0^h \left( f(x)+(t-x)f'(x)+\dfrac{1}{2}(t-x)^2f''(x)+\cdots \right)\mbox{d}t \\[0.8em] & = & hf(x)+h\left(\dfrac{h}{2}-x\right)f'(x)+O(h^3) \end{array} \] de modo que \[ \tag{3.12} \mbox{Viés}\big(\widehat{f}(x)\big) = \dfrac{p_0}{h}-f(x)=\left(\dfrac{h}{2}-x\right)f'(x)+O(h^2)\cdot \] Consulte a Figura 3.3. Para referência futura, observe que (3.12) implica que o viés é de ordem superior \(O(h^2)\) no centro de um intervalo, quando \(x = h/2\).

Figura 3.3: Viés do estimador de histograma em um intervalo típico.

Usando o teorema do valor médio generalizado (GMVT), o termo principal do viés quadrático integrado (ISB) para este intervalo é \[ \tag{3.13} \int_{B_0} \left(\dfrac{h}{2}-x\right)^2 f'(x)^2\mbox{d}x = f'(\eta_0)^2\int_0^h \left(\dfrac{h}{2}-x\right)^2\mbox{d}x = \dfrac{h^3}{12}f'(\eta_0)^2, \] para algum \(\eta_0\in B_0\).

Este resultado se generaliza para outros compartimentos para alguma coleção de pontos \(\eta_k\in B_k\). Portanto, o ISB total é \[ \tag{3.14} \mbox{ISB} = \dfrac{h^2}{12}\sum_{k=-\infty}^\infty f'(\eta_k)^2\ h = \dfrac{h^2}{12}\int_{-\infty}^\infty f'(x)^2\mbox{d}x+o(h^2), \] o que decorre da convergência Riemanniana padrão de somas para integrais.

Uma nota sobre as hipóteses: se a variação total de \((f')^2(\cdot)\) for finita, então o termo residual em (3.14) torna-se \(O(h^3)\) (ver Exercício 3.3). Assumindo a existência de uma segunda derivada absolutamente contínua, o termo de erro é \(O(h^4)\). Há pouca diferença prática entre essas observações.

Em vez de insistir em acompanhar explicitamente a ordem do termo residual, a seguinte notação será adotada. Os termos principais no ISB serão referidos como o ISB assintótico (AISB). Assim, \(\mbox{ISB} = \mbox{AISB} + o(h^2)\) e \[ \mbox{AISB} = \dfrac{1}{12}h^2 \int_{-\infty}^\infty f'(x)^2\mbox{d}x=\dfrac{1}{12}h^2 R(f')\cdot \] De forma semelhante, a variância integrada assintótica (AIV) e o erro quadrático médio integrado assintótico (AMISE) referem-se aos termos principais nas aproximações da IV e do MISE, respectivamente. O teorema a seguir, que resume as equações (3.11) e (3.14).


Teorema 3.3 Suponha que \(f\) tenha uma derivada absolutamente contínua e uma primeira derivada de quadrado integrável. Então, o MISE assintótico é \[ \mbox{AMISE}(h) = \dfrac{1}{nh}+\dfrac{1}{12}h^2 R(f'), \] portanto \[ h^* = \Big(6/R(f') \Big)^{1/3} n^{-1/3} \] e \[ \tag{3.15} AMISE^* = (3/4)^{2/3} \big(R(f')\big)^{1/3} n^{-2/3}\cdot \]


Demonstração Scott (1979) e Freedman and Diaconis (1981).

Assim, a largura do intervalo assintoticamente ótima depende da densidade desconhecida apenas por meio da rugosidade de sua primeira derivada. Esse resultado se mantém independentemente da escolha da origem do intervalo, que deve ter um papel secundário no MISE em comparação com a largura do intervalo.

A largura ótima do intervalo diminui a uma taxa relativamente lenta. O erro ótimo correspondente, \[ \mbox{AMISE}^* = \mbox{AMISE}(h^*), \] diminui na mesma taxa que o limite no Corolário 3.2, longe da taxa desejável \(O(n^{-1})\).


3.2.3 Regra de referência da densidade Normal


Será conveniente fixar alguns desses resultados para o caso especial de dados normais. Se \[ f=N(\mu,\sigma^2), \] então, \(R(f')=1/\big(4\sqrt{\pi}\sigma^3 \big)\) (veja Exercício 3.2). Portanto, a partir do Teorema 3.3, \[ \tag{3.16} h^*=\big(24\sqrt{\pi} \sigma^3/n\big)^{1/3} \approx 3.5 \, \sigma \, n^{-1/3}\cdot \] Scott (1979) propôs usar a densidade normal como densidade de referência para construir histogramas a partir de dados amostrados, usando o desvio padrão amostral \(\widehat{\sigma}\) em (3.16) para obter o \[ \tag{3.17} \mbox{Regra de referência para largura normal do intervalo:}: \quad \widehat{h}= 3.5\, \widehat{\sigma} \, n^{-1/3}\cdot \]

Como \(\widehat{\sigma}\to\sigma\) ocorre mais rapidamente que \(O(n^{-1/3})\), essa regra é muito estável. Freedman and Diaconis (1981) propuseram uma regra mais robusta, substituindo o parâmetro de escala desconhecido \(\sigma\) por um múltiplo do intervalo interquartil (IQR): \[ \tag{3.18} \widehat{h}=2 (\mbox{IQR})n^{-1/3}\cdot \] Se os dados forem de fato normais, a regra de Freedman-Diaconis corresponde a cerca de 77% da regra de Scott, visto que o \(\mbox{IQR} = 1.348 \sigma\) neste caso.


3.2.3.1 Comparação das regras de largura de banda


Como a regra ótima se compara à regra de Sturges no caso em que a situação é mais favorável a esta última? Na Tabela 3.1, é apresentado o número de intervalos sugeridos pelas três regras, assumindo que os dados seguem uma distribuição normal e que o intervalo da amostra é \((-3,3)\).

Tabela 3.1: Comparação do número de intervalos de três regras de referência normais.

Talvez a observação mais notável que se possa obter desta tabela seja a grande concordância entre as regras para amostras entre 50 e 500 pontos. Para amostras maiores, a regra de Sturges fornece um número insuficiente de intervalos, resultando em uma estimativa de histograma excessivamente suavizada que desperdiça grande parte da informação contida nos dados. Como mostrado na Seção 3.3.1, praticamente qualquer outra densidade amostral não normal exigirá ainda mais intervalos.

Este resultado é consequência do fato de que o caso normal está muito próximo de um limite inferior teórico para o número de intervalos ou, equivalentemente, de um limite superior para a largura do intervalo. A regra de Freedman-Diaconis apresenta 35% mais intervalos do que a regra de Scott e o histograma será mais irregular.

Figura 3.4: AMISE versus largura de banda para uma densidade \(\mbox{Beta}(5,5)\). As melhores larguras de banda e as de Sturges são indicadas por pontos nas curvas. As larguras de banda de referência de Freedman-Diaconis e Scott são mostradas como pontos semicirculares e triangulares, respectivamente, ao longo do eixo \(x\).

Como outro exemplo, considere a densidade beta, \(\mbox{Beta}(5,5) = 630x^4 (1-x)^4\), que tem suporte finito em \((0,1)\) e é semelhante a uma densidade normal. As curvas assintóticas de MISE versus largura de banda são mostradas na Figura 3.4 para vários tamanhos de amostra. Novamente, a regra de Sturges leva a um suavização excessiva severa à medida que o tamanho da amostra aumenta.

Também são mostradas ao longo do eixo as regras de referência de Scott e Freedman-Diaconis. A regra de Freedman-Diaconis (robusta) é calibrada para ser mais agressiva do que uma regra de referência puramente normal, sendo 18% mais estreita do que a regra de referência de Scott para esta densidade.

O IQR e o desvio padrão são 0.16 e \(1/\sqrt{44} = 0.151\), de modo que as constantes de largura de banda nas Equações (3.17) e (3.18) são \(2 \, \mbox{IQR} = 0.432\) e \(3.5\, \sigma = 0.528\), respectivamente, uma diferença de 18%.


3.2.3.2 Ajustes para assimetria e curtose


Scott (1979) também propôs o uso das distribuições lognormal e \(t\)-Student como densidades de referência para modificar a regra normal quando os dados são assimétricos ou apresentam caudas pesadas. Suponha que \(Y\) seja uma variável aleatória não normal com densidade \(g(y)\) e momentos \(\sigma_Y^2\), \(\beta_1\) e \(\beta_2\) a assimetria e a curtose padronizadas, respectivamente.

Considere o uso da regra normal com \(\sigma = \sigma_Y\) em comparação com a largura de intervalo ótima no Teorema 3.3. A razão entre essas duas larguras de intervalo é dada por: \[ \tag{3.19} \dfrac{h_Y^*}{h_N}=\left(\dfrac{R\Big(\phi'\big(y\, | \, 0,\sigma_Y^2\big)\Big)}{R\big(g'(y)\big)} \right)^{1/3}; \] se essa razão for muito diferente de 1, então a regra de referência normal deve ser modificada por essa razão.

Por exemplo, seja \(g(y)\) a densidade lognormal da variável aleatória \(Y = \exp(X)\), onde \(X\sim N(0,\sigma^2)\). Então, tratando \(\sigma^2\) como um parâmetro, temos: \[ \sigma_Y^2=e^{\sigma^2}\big(e^{\sigma^2}-1\big); \quad \beta_1=\big(e^{\sigma^2}-1\big)^{1/2}\big(e^{\sigma^2}+2\big); \quad R(g')=\dfrac{(\sigma^2+2)e^{9\,\sigma^2/4}}{8\sqrt{\pi}\, \sigma^3}\cdot \]

Como \(R(\phi'(y)) = 1/\big(4\, \sqrt{\pi}\, \sigma_Y^3\big)\), a regra normal deve ser multiplicada pelo fator em (3.19) dado por \[ \tag{3.20} \mbox{Fator de assimetria}\big(\beta_2(\sigma))\big)=\dfrac{2^{1/3}\sigma}{e^{4 \, \sigma^2/4}\big(\sigma^2+2\big)^{1/3} \big(e^{\sigma^2}-1 \big)^{1/2} }\cdot \]

Este fator é apresentado na Figura 3.5. Claramente, qualquer assimetria requer larguras de intervalo menores do que as fornecidas pela regra de referência normal.

Figura 3.5: Fatores que modificam as regras normais de largura do intervalo de referência para o histograma e o polígono de frequência em função da assimetria padronizada e da curtose excessiva.

Um cálculo semelhante pode ser realizado para a curtose, assumindo que \(Y\sim t_\nu\), que é a distribuição \(t\) com \(\nu\) graus de liberdade. Seja o excesso de curtose denotado por \(\widetilde{\beta}_2 = \beta_2 - 3\). Tratando \(\nu\) como um parâmetro, obtemos: \[ \sigma_Y^2=\dfrac{\nu}{\nu-2}; \quad \widetilde{\beta}_2=\dfrac{6}{\nu-4}; \quad R(g')=\dfrac{2\, \Gamma(\nu+3/2)\Gamma^2\big((\nu+3)/2\big)}{\sqrt{\pi}\, \nu^{3/2} \Gamma(\nu+3)\Gamma^2(\nu/2)}\cdot \] Substituindo em (3.19), obtemos \[ \tag{3.21} \mbox{Fator de kurtosis}\Big(\widetilde{\beta}_2(\nu) \Big)=\dfrac{\sqrt{\nu-2}}{2}\left(\dfrac{\Gamma(\nu+3)\Gamma^2(\nu/2)}{\Gamma(\nu+3/2)\Gamma^2\big((\nu+3)/2\big)} \right)^{1/3}\cdot \] Veja a Figura 3.5. A modificação na largura do intervalo não é tão grande quanto para alta assimetria.


3.2.4 Tamanhos de amostra equivalentes


Em comparação com estimadores paramétricos, os histogramas têm a vantagem de uma consistência bastante geral, sem o que Fisher chamou de “problema de especificação”. Qual é a penalidade incorrida se um histograma ótimo for usado em vez do estimador paramétrico ótimo correspondente?

Considere novamente o exemplo da Seção 3.2.3. Do Teorema 3.3, segue-se que se \(f = N(0,1)\), então o \(\mbox{AMISE}\) ótimo do histograma é \[ \mbox{AMISE}^*=\Big(9/\big(64\sqrt{\pi}\big) \Big)^{1/3}n^{-2/3}\approx 0.4297 \, n^{-2/3}\cdot \] Suponha que uma amostra de tamanho 100 esteja disponível para o estimador paramétrico \(N(\overline{x},\sigma^2)\), para o qual \(\mbox{AMISE}\) = 0.002468 da equação (2.9). A Tabela 3.2 apresenta os tamanhos de amostra equivalentes para o histograma e vários estimadores paramétricos. Os estimadores paramétricos são claramente superiores, principalmente para erros menores.

Tabela 3.2: Tamanhos de amostra equivalentes para vários estimadores da densidade normal.

Por outro lado, se a densidade verdadeira for apenas aproximadamente normal, o \(\mbox{MISE}\) paramétrico nunca será menor que o nível de viés quadrático integrado verdadeiro. \[ \mbox{ISB}=\int \Big( \phi\big( x\, | \, \mu_f,\sigma_f^2\big) - f(x)\Big)^2 \mbox{d}x, \] onde \(\mu_f\) e \(\sigma_f^2\) são os momentos reais da densidade desconhecida \(f(x)\).

O estimador de máxima verossimilhança, assumindo normalidade, corresponderá assintoticamente a esses momentos. Assim, embora a variância dos parâmetros tenda a zero assintoticamente, o \(\mbox{ISB}\) (Índice de Estabilidade de Coeficientes) permanecerá. Testar a qualidade do ajuste de um modelo paramétrico é necessário e pode ser especialmente difícil no caso multivariável; distinguir entre certas famílias paramétricas pode exigir amostras grandes.


3.2.5 Sensibilidade do MISE à largura do intervalo


Qual a importância de se escolher o parâmetro de suavização de forma otimizada? Se praticamente qualquer largura de intervalo funcionar, então será necessário menos esforço para resolver o problema de calibração.


3.2.5.1 Caso assintótico


Na prática, dados dados reais de uma densidade desconhecida, o parâmetro de suavização escolhido não será \(h^∗\), mas sim da forma \(h = ch^∗\).

Em média, se \(c\ll 1\), o histograma terá alta variância e será “muito irregular”; se \(c\gg 1\), a estimativa terá alto viés ou erro sistemático e será “muito suave”. Quão sensível é o MISE a desvios locais de \(c\) em relação a 1?

O MISE assintótico do histograma no Teorema 3.3 tem a forma \[ \tag{3.22} \mbox{AMISE}(h)=\dfrac{a}{nh}+\dfrac{b}{2}h^2, \] onde \(a\) e \(b\) são constantes positivas.

Em vez de considerar apenas este caso específico, será analisada uma forma mais geral do AMISE: \[ \tag{3.23} \mbox{AMISE}(h)=\dfrac{a}{(d+2r)\, nh^{d+2r}}+\dfrac{b}{2p}h^{2p}, \] onde \((d,p,r)\) são inteiros positivos e \(a\) e \(b\) são constantes positivas que dependem do estimador de densidade e da função de densidade desconhecida.

Em geral, a tripla \((d,p,r)\) refere-se, respectivamente:

    1. à dimensão dos dados,

    2. à “ordem” do viés do estimador e

    3. à ordem da derivada que está sendo estimada.

Comparando as equações (3.22) e (3.23), \((d,p,r) = (1,1,0)\) para o caso do histograma.

Da minimização da equação (3.23), segue-se que \[ h^*=\big(a/b\, n\big)^{1/d+2r+2p} \] e \[ \tag{3.24} \mbox{AMISE}^* = \mbox{AMISE}(h^*)=\left(\dfrac{d+2p+2r}{2p(d+2r)}\right)\left(\dfrac{a^{2p}b^{d+2r}}{n^{2p}}\right)^{1/(d+2p+2r)}\cdot \] No Exercício 3.5, demonstra-se que a porção da variância do \(\mbox{AMISE}^*\) do histograma é o dobro da porção do viés ao quadrado.

Finalmente, é pode-se verificar (ver Exercício 3.6) que \[ \tag{3.25} \dfrac{\mbox{AMISE}(ch^*)}{\mbox{AMISE}(h^*)}=\dfrac{2p+(d+2r)c^{d+2p+2r}}{(d+2p+2r)c^{d+r}}, \] que, no caso do historama é igual a \((2+c^2)/3c\).

Uma consequência dessa expressão é que os desvios de \(h\) em relação a \(h^∗\) devem ser medidos de forma multiplicativa e não aditiva. Isso também pode ser observado em uma análise dimensional de \(h^∗\) no Teorema 3.3. Na experiência do autor, uma mudança de 10-15% em \(h\) produz uma mudança pequena, porém perceptível, no histograma.

No entanto, na Tabela 3.3, fica claro que o histograma, que corresponde ao caso \(p = 1\), é bastante insensível à escolha da largura do intervalo dentro de 33% do valor ótimo. Observe que o critério \(L_2\) é menos afetado por erros de alta variância do que por erros de alto viés, por exemplo, \(c = 1/2\) vs. \(c = 2\).

Tabela 3.3: Sensibilidade do AMISE ao erro na escolha da largura do intervalo \(h = ch^∗\).

Olhando para o futuro, valores maiores de \(p\) correspondem a métodos de “ordem superior” que apresentam taxas de MISE mais rápidas. No entanto, a sensibilidade à escolha dos parâmetros de suavização é muito maior. O aumento da sensibilidade também é evidente com o aumento da dimensão \(d\) e da ordem da derivada \(r\) (ver Exercício 3.7).


3.2.5.2 Simulações com amostras grandes e pequenas


A análise anterior focou apenas no comportamento médio. O ISE real para uma amostra individual é considerado nesta seção. A Figura 3.6 exibe três histogramas de 1000 observações normais. Quando \(h = 2h^*\), a estimativa apresenta um viés substancial devido à largura dos intervalos. Em consultas informais, a maioria dos estatísticos prefere o histograma com \(h = \frac{1}{2}h^*\), mesmo que ele contenha vários pequenos picos espúrios. No entanto, é fácil suavizar visualmente o ruído local nessa estimativa. O inverso não é verdadeiro, já que não é possível visualizar os detalhes perdidos no histograma com largura de intervalo \(h = 2h^*\).

Figura 3.6: Três histogramas de 1000 observações normais com larguras de banda ou largura de intervalo \(h = \big(\frac{1}{2}h^*,h^*,2h^*\big)\).

Generalizar demais a partir de experiências com amostras pequenas é uma tentação que deve ser evitada. A Figura 3.7 exibe quatro histogramas de uma amostra de um milhão de pontos normais. O \(\mbox{ISE}\) exato é fornecido em cada quadro.

Figura 3.7: \(\mbox{ISE}\) exato de quatro histogramas de um milhão de pontos normais com larguras de intervalo \(h=(h^*/2,h^*,2h^*,4h^*)\), onde \(h^*=0.035\).

Examine as mudanças no \(\mbox{ISE}\). Localmente, o histograma permanece muito sensível a mudanças em \(h\), mesmo quando \(n\) é grande. Isso pode parecer contraintuitivo inicialmente, especialmente considerando as condições muito gerais necessárias para a consistência. Mas amostras grandes não garantem automaticamente boas estimativas.

O histograma com \(h = \frac{1}{2}h^∗\) contém muitas modas locais. O comportamento das modas amostrais em um histograma será reexaminado na Seção 3.5.2. Mesmo o histograma “ótimo” com \(h = h^∗\) parece ruidoso (localmente) em comparação com o histograma com \(h = 2h^∗\). A observação de que estimativas suavizadas otimamente parecem estar à beira de uma forma de instabilidade parece ser uma ocorrência comum.


3.2.6 MISE exato versus MISE assintótico


Para uma malha geral com espaçamento desigual, onde a largura do intervao \(B_k\) é \(h_k\), podemos mostrar a partir das equações (3.3) e (3.4) que \[ \tag{3.26} \mbox{IV}=\dfrac{1}{n}\sum_k \dfrac{p_k(1-p_k)}{h_k} \qquad \mbox{e} \qquad \mbox{ISB}=R(f)-\sum_k \dfrac{p_k^2}{h_k}, \] exatamente (veja o exercício 3.8).

Agora, \(\sum_k p_k=1\); portanto, para o caso especial de uma malha igualmente espaçada, \(h_k = h\), \[ \mbox{MISE}(h,t_0,n)=\dfrac{1}{nh}-\dfrac{n+1}{nh}\sum_k p_k^2+R(f)\cdot \] Esta expressão pode ser usada para testar a precisão do \(\mbox{AMISE}\) para amostras pequenas, para verificar a sensibilidade do \(\mbox{MISE}\) à especificação incorreta da largura do intervalo e para avaliar o efeito da escolha da origem da malha \(t_0\).


3.2.6.1 Densidade Normal


Na Figura 3.8, o \(\mbox{MISE}\) exato e assintótico de um histograma de largura de banda fixos com \(t_0 = 0\) são mostrados graficamente em escala log-log para dados normais padrão com vários tamanhos de amostra. Claramente, o \(\mbox{AMISE}\) aproxima adequadamente o \(\mbox{MISE}\) verdadeiro, mesmo com amostras pequenas.

O erro de aproximação (lacuna) diminui rapidamente à medida que \(n\) aumenta. De particular interesse é o fato de que as larguras das bandas que minimizam o \(\mbox{MISE}\) e o \(\mbox{AMISE}\) são praticamente indistinguíveis para todos os valores de \(n\). Que tanto \(h\) quanto o \(\mbox{MISE}\) diminuem à medida que \(n\) aumenta é facilmente perceptível. Além disso, a lacuna é explicada quase inteiramente pelo termo \(-R(f)/n\) da equação (3.11).

Figura 3.8: \(\mbox{AMISE}\) e \(\mbox{MISE}\) exato para a densidade \(N(0,1)\).

Na Figura 3.9, é apresentada a decomposição do \(\mbox{MISE}\) em \(\mbox{IV}\) e \(\mbox{ISB}\). As linhas de \(\mbox{IV}\) e \(\mbox{ISB}\) na escala log-log são quase retas, com inclinações conforme previsto pela teoria assintótica no Teorema 3.3.

Figura 3.9: Decomposição do \(\mbox{MISE}\) para a densidade \(N(0,1)\) do viés quadrático integrado/variância.

Talvez o mais surpreendente seja a linha única para o \(\mbox{ISB}\). Essa característica é facilmente explicada pelo fato de que o \(\mbox{ISB}\), conforme dado na equação (3.26), é função apenas da partição e não do tamanho da amostra. Observe que, na escolha ótima da largura do intervalo, a contribuição do \(\mbox{IV}\) excede a do \(\mbox{ISB}\) em uma proporção aproximada de 2:1, como mostrado no exercício 3.5. O critério \(L_2\) é mais sensível quando \(h > h^*\).

A Figura 3.10 exibe uma imagem da curva \(\mbox{MISE}\) em uma faixa muito mais ampla de larguras de intervalo quando \(n = 100\). O \(\mbox{MISE}(h)\) é ilimitado quando \(h\to 0\), mas \[ \lim_{h\to\infty} \mbox{MISE}(h) =\lim_{h\to\infty} \mbox{ISB}(h)=R(f), \] como é claramente evidente.

Figura 3.10: Curva \(\mbox{MISE}\) completa para a densidade \(N(0,1)\) quando \(n = 100\). A relação de sensibilidade assintótica é válida em uma ampla gama de larguras de intervalo.

A teoria assintótica prevê a maneira como o \(\mbox{MISE}\) varia com \(h\) (ver equação (3.25)). A linha tracejada na Figura 3.10 é um gráfico de \[ \left(h,\dfrac{c^3+2}{3\, c}\mbox{MISE}^*\right), \] onde \(c=h/h^*\) e \(\mbox{MISE}^*=\mbox{MISE}(h^*)\). Essa aproximação é precisa muito além da vizinhança imediata de \(h^∗\).


3.2.6.2 Densidade Lognormal


A Figura 3.11 repete o cálculo indicado na Figura 3.8, mas com uma densidade Lognormal e \(t_0 = 0\). Embora a \(\mbox{AMISE}\) mantenha seu formato parabólico, as curvas \(\mbox{MISE}\) reais são bastante complexas para \(n < 1000\).

Figura 3.11: \(\mbox{AMISE}\) e \(\mbox{MISE}\) exato para a densidade Lognormal.

Quando \(n = 25\), a teoria assintótica prevê \(h^∗ = 0.49\), enquanto \(h = 1.28\) é, na verdade, o valor ótimo. Além disso, \(\mbox{AMISE}^∗ = 0.1218\), que é muito maior que o real \(\mbox{MISE}^∗=0.0534\). Ademais, quando \(n = 200\), a curva MISE apresenta 2 (!) mínimos locais em \(h = 0.245\) e \(h = 0.945\). Para \(n > 1000\), a teoria assintótica parece ser adequada.

O que pode explicar essas observações? Primeiro, a densidade Lognormal, embora infinitamente diferenciável, é estatisticamente muito irregular, particularmente perto da origem. De fato, 90% da irregularidade, \(R(f') = 3\, e^{9/4}/\big(8\sqrt{\pi}\big)\), vem do pequeno intervalo [0,0.27], mesmo que a moda esteja localizada em \(x = 0.368\) e o percentil 99 seja \(x = 10.2\).

Portanto, assintoticamente, a largura ideal do intervalo deve ser relativamente pequena para acompanhar a densidade que muda rapidamente perto da origem. Para amostras muito pequenas, não há dados suficientes para acompanhar com precisão o aumento da densidade Lognormal à esquerda da moda; daí a largura de banda ideal do \(\mbox{MISE}\) \(h > 1\), que é muito maior do que a prevista pelo \(\mbox{AMISE}\). Para \(n > 1000\), a largura ideal do intervalo satisfaz \(h < 0.15\) e o aumento pode ser acompanhado. Para \(n\approx 200\), as duas larguras de intervalo localmente ótimas indicam uma situação em que a subida pode ser igualmente bem aproximada por um intervalo estreito ou largo, ou por quase qualquer largura de intervalo intermediária, uma vez que a curva MISE é relativamente plana.

A Figura 3.11 ilustra uma observação prática e generalizável importante sobre métodos não paramétricos. Os tamanhos de amostra podem ser agrupados em três categorias: inadequados, de transição e suficientes.

Para dados Lognormais, os tamanhos de amostra inadequados correspondem a \(100 < n < 400\) com histogramas igualmente espaçados e os de transição a \(400 < n < 1000\). Para outras densidades, a região de transição pode se estender a amostras ainda maiores, embora a densidade Lognormal seja razoavelmente extrema. Infelizmente, na prática, com um conjunto de dados específico, não é possível ter certeza em qual categoria \(n\) se enquadra.

No entanto, como pode ser facilmente observado, características mais estreitas que \(h\) não podem ser detectadas de forma confiável abaixo de certos tamanhos de amostra. Donoho (1988) resume de forma concisa a situação na estimação não paramétrica como inferência unilateral — somente tamanhos de amostra maiores podem responder definitivamente se existe uma estrutura menor na densidade desconhecida.


3.2.7 Influência da localização da borda do intervalo no MISE


3.2.7.1 Caso geral


A análise assintótica do MISE no Teorema 3.3 indica que a escolha da origem da malha de intervalos \(t_0\) tem um efeito de ordem inferior, desde que as condições do teorema sejam satisfeitas. Para dados \(N(0,1)\) com suavização ótima, o MISE exato é minimizado se o ponto \(x = 0\) estiver no meio de um intervalo, em vez de na borda. A diferença no MISE ao usar essas duas localizações da borda dos intervalos é minúscula, sendo de 1.09% quando \(n = 25\) e menor que \(10^{-5}\) quando \(n > 100\).

Com dados Lognormais, o efeito de borda também é menor que \(10^{-5}\) quando \(n > 400\). Mas quando \(n = 25\) e \(h = 1.28\), o MISE varia entre (0.042, 0.149) — o melhor para a malha \[ (−1.20, 0.08, 1.36, 2.64,\cdots) \] e o pior para a malha \[ (−0.65, 0.63, 1.91, 3.19,\cdots)\cdot \]

Claramente, a escolha da borda do intervalo é desprezível para tamanhos de amostra suficientemente grandes, mas pode ser significativa para tamanhos de amostra inadequados. Graficamente, a escolha de \(t_0\) para \(h\) fixo fornece uma gama de histogramas possíveis com aparências subjetivas bastante diferentes (veja o Capítulo 5 e a Figura 5.1).


3.2.7.2 Descontinuidades de fronteira na densidade


As propriedades de aproximação de um histograma não são afetadas por um salto simples na densidade, se o salto ocorrer na fronteira de um intervalo do histograma. Isso decorre de uma análise cuidadosa da derivação do \(\mbox{MISE}\) na Seção 3.2.2. Descontinuidades podem afetar negativamente todos os estimadores de densidade. O efeito adverso pode ser demonstrado para o histograma por meio de um exemplo.


Exemplo 3.1: Histograma na desidade Exponencial.

Considere \(f(x) = e^{-x}\), \(x\geq 0\) e a malha \(t_k = kh\) para inteiro \(k\geq 0\). Então, o Teorema 3.3 é válido no intervalo \((0,\infty)\), no qual \(R(f') = 1/2\). Portanto, \[ h^*=\big(12/n\big)^{1/3} \qquad \mbox{e} \qquad \mbox{AMISE}^* = 0.6552\, n^{-2/3}\cdot \]

Suponha que a descontinuidade em zero não fosse conhecida a priori e a malha \(t_k = \big(k-\frac{1}{2}\big)h\) fosse escolhida. Então, a atenção se concentra no viés no intervalo \(B_0 = [-h/2, h/2)\) onde a densidade é descontínua (veja a Figura 3.12).

Figura 3.12: Ilustração do problema de descontinuidade de fronteira em compartimentos.

Observe que a massa de probabilidade no compartimento \(B_0=(-h/2,h/2]\) é \[ p_0=\int_0^{h/2} e^{-x}\mbox{d}x=1-e^{-h/2}\cdot \] Agora, \(\mbox{E}\big(\widehat{f}(x) \big)=p_0/h\); portanto, \[ \tag{3.27} \int_{-h/2}^{h/2} \mbox{Viés}^2(x)\mbox{d}x =\int_{-h/2}^0 \left(\dfrac{p_0}{h}-0 \right)^2\mbox{d}x + \int_0^{h/2} \left(\dfrac{p_0}{h}-e^{-x} \right)^2\mbox{d}x, \] que é igual a \[ \dfrac{1-e^{-h}}{2}+\dfrac{2\, e^{-h/2}-e^{-h}-1 }{h}\approx \dfrac{h}{4}-\dfrac{h^2}{8}+\cdots \, \cdot \] Assim, ao longo do intervalo \((-h/2,\infty)\), \[ \tag{3.28} \mbox{ISB}=\dfrac{h}{4}-\dfrac{h^2}{8}+O(h^3)+\dfrac{h^2}{12}\int_{h/2}^{\infty} f'(x)^2\mbox{d}x+o(h^2)\cdot \] O pior resultado foi alcançado — o ISB é inteiramente dominado pela contribuição do intervalo que contém a descontinuidade em \(x = 0\). A integral da variância assintótica permanece inalterada, de modo que \[ \mbox{AMISE}(h)=\dfrac{1}{nh}+\dfrac{h}{4} \] o que implica \[ \tag{3.29} h^*=\dfrac{2}{\sqrt{n}} \qquad \mbox{e} \qquad \mbox{AMISE}^*=\dfrac{1}{\sqrt{n}}, \] o que é significativamente pior do que \(O(n^{-2/3})\).

Na verdade, a taxa é tão lenta quanto para dados bivariados, como será mostrado na Seção 3.4. A Tabela 3.4 ilustra os custos. O histograma tenta acomodar a descontinuidade escolhendo uma largura de intervalo menor. Compare essa situação com o comportamento muito irregular da densidade Lognormal perto da origem com amostras pequenas.

Tabela 3.4: Impacto potencial no \(\mbox{AMISE}\) da falta de conhecimento das descontinuidades de fronteira.


3.2.8 Malhas de histograma otimamente adaptativas


Um histograma com largura de intervalo fixa, calibrado de forma otimizada com dados reais, frequentemente apresenta irregularidades nas caudas devido à escassez de dados. Esse fenômeno é uma das várias razões para se considerar histogramas com malhas adaptativas.

Avaliar a redução do \(\mbox{MISE}\) com malhas adaptativas é a tarefa desta seção. Criar malhas adaptativas é um exercício familiar para qualquer pessoa que já tenha realizado um teste de aderência \(\chi^2\): para atender à recomendação de que a contagem esperada em cada célula exceda 5, as células nas caudas geralmente são combinadas ou, alternativamente, as células são formadas de modo que cada uma contenha exatamente o mesmo número de pontos (ver Seção 3.2.8.4).


3.2.8.1 Limites para a melhoria do MISE em histogramas adaptativos


Um limite inferior, assintoticamente correto, para a redução do MISE de um histograma otimamente adaptativo decorre das equações (3.6) e (3.14). Vamos construir uma aproximação para o MSE pontual do histograma assintoticamente adaptativo (AAMSE): \[ \tag{3.30} \mbox{AAMSE}(x)\approx \dfrac{f(x)}{nh}+\dfrac{h^2}{12} f'(x)^2\cdot \] Minimizando o \(\mbox{AAMSE}(x)\) para cada ponto \(x\). Portanto, \[ h^*(x)=\left(\dfrac{6\, f(x)}{nf'(x)^2} \right)^{1/3} \] e, então, \[ \mbox{AAMSE}^*(x)=\left(\dfrac{3\, f(x)f'(x)}{4n} \right)^{2/3}\cdot \] Integrando \(\mbox{AAMSE}^∗(x)\) em relação a \(x\), obtemos o seguinte resultado.


Teorema 3.4: Assintoticamente, para um histograma otimamente adaptativo, \[ \tag{3.31} \displaystyle \mbox{AAMISE}^*=(3/4)^{2/3} \left(\int_{-\infty}^\infty \Big(f'(x)f(x) \Big)^{2/3}\mbox{d}x \right) \, n^{-2/3}\cdot \]


Demonstração Este resultado foi discutido por Terrell and Scott (1983), Terrell and Scott (1992) e Kogure (1987).

Comparando as equações (3.15) e (3.31), a melhoria da malha adaptativa é garantida se \[ \int \Big(f'(x)f(x) \Big)^{2/3}\mbox{d}x \leq \left(\int f'(x)^2\mbox{d}x \right)^{1/3} \] ou, equivalentemente, se \[ \tag{3.32} \mbox{E}\left(\left(\dfrac{f'(X)^2}{f(X)} \right)^{1/3}\right)\leq \left(\mbox{E} \left(\dfrac{f'(X)^2}{f(X)}\right)\right)^{1/3}; \] mas esta última desigualdade decorre da desigualdade de Jensen (a versão da função côncava) (ver exercício 3.11).

A Tabela 3.5 mostra cálculos numéricos dessa razão para algumas densidades simples. À primeira vista, os ganhos com malhas adaptativas são surpreendentemente pequenos.

Tabela 3.5: \(\mbox{AMISE}\) reduzido usando uma malha de histograma otimizada e adaptativa.


3.2.8.2 Algumas malhas ótimas


Como a expressão para o MISE exato para malhas arbitrárias está disponível na equação (3.26), a malha adaptativa ótima pode ser encontrada numericamente. Certas características nessas malhas ótimas podem ser previstas, veja a discussão após a equação (3.8). Comparadas à largura de banda fixa ótima, as classes serão relativamente mais largas não apenas nas caudas, mas também perto de modas onde \(f'(x)\) é pequeno. Em regiões onde a densidade está mudando rapidamente, as classes serão mais estreitas.

Esse padrão de “acordeão” é aparente no exemplo a seguir com uma densidade \(Beta(5,5)\) transformada para o intervalo \((-1, 1)\), ou seja, \[ f(x) = \dfrac{315}{256} (1-x^2)^4_+, \] (veja a Figura 3.13).

Como a função de distribuição acumulada (FDA) dessa densidade é um polinômio, ela é mais adequada para otimização numérica do que a FDA normal, em comparação com o MISE adaptativo exato dado na equação (3.26). Claramente, a malha ótima é simétrica em relação a 0 neste exemplo. Mas a malha ótima não inclui 0 como um nó (nó da malha). Forçar 0 a ser um nó aumenta o MISE, por exemplo, em 4.8% quando \(n = 10.000\); o melhor MISE para uma malha com espaçamento uniforme \(h = 0.0515\) é 8.4% maior. Esses resultados são muito semelhantes aos da densidade normal.

Figura 3.13: Representação das malhas adaptativas ótimas para a densidade \(Beta(5,5)\) transformada, que é igual a \(315/256(1-x^2)^4_+\). A malha adaptativa ótima também é indicada por marcas de escala acima de cada gráfico.


3.2.8.3 Espaço nulo de densidades adaptativas


Existe uma classe de funções de densidade para as quais não há melhoria assintótica com um procedimento de histograma adaptativo. Essa classe pode ser chamada de espaço nulo de densidades adaptativas.

Analisando a equação (3.32), qualquer densidade que resulte em igualdade na desigualdade de Jensen está no espaço nulo. Isso ocorre quando o argumento entre colchetes é constante; ou seja, \[ \dfrac{f'(x)^2}{f(x)}=c, \] o que implica \[ \tag{3.33} f(x)=\dfrac{1}{4}c \, (x-a)^2, \] onde “a” é uma constante arbitrária.

As densidades \(f_1(x) = 3(1-x)^2 \pmb{I}_{[0,1]}(x)\) e \(f_2(x) = \frac{3}{2}(1-|x|^2)\pmb{I}_{[-1,1]}(x)\) também satisfazem a equação diferencial por partes. De fato, qualquer densidade obtida pela junção de densidades da forma (3.33) de maneira contínua está no espaço nulo, veja a Figura 3.14.

Figura 3.14: Exemplos de densidades quadráticas por partes para as quais a malha adaptativa assintoticamente ótima é, na verdade, igualmente espaçada.

Como antes, a melhor malha adaptativa pode ser encontrada por otimização numérica. No entanto, com \(n = 100\) e a densidade \(f_1(x)\) definida anteriormente, os nós da malha ótima calculada diferem em menos de 2% da malha igualmente espaçada com \(h = 1/6\). E o \(\mbox{MISE}\) da malha ótima calculada é apenas 0.06% melhor do que o da malha igualmente espaçada.


3.2.8.4 Malhas percentis


Na prática, encontrar a malha adaptativa descrita na seção anterior é difícil. Assim, alguns autores propuseram malhas adaptativas (não ótimas) ou histogramas adaptativos com contagens iguais de pontos em cada classe, que são mais fáceis de implementar. Uma das malhas mais intuitivamente atraentes possui um número igual (ou fração igual) de pontos em cada classe.

Essa ideia pode ser modelada por uma malha percentil, que é não estocástica, com \(m\) compartimentos, cada um contendo uma fração idêntica de massa de probabilidade igual a \(1/m\). A malha terá \(m+1\) nós em \[ \tag{3.34} t_k = F^{-1}_X\Big(\dfrac{k}{m} \Big), \qquad k=0,1,\cdots,m\cdot \]

Antes de calcular o MISE exato para tais malhas, é instrutivo apresentar uma derivação heurística. Assintoticamente, para qualquer \(x\in B_i\), \(h_i f(x)\approx p_i\), a massa de probabilidade no \(i\)-ésimo compartimento de largura \(h_i\). A equação (3.30) sugere que \[ \mbox{AAMSE}(x)\approx \dfrac{f(x)}{nh_i}+\dfrac{h^2_i f'(x)^2}{12}\cdot \] Se cada compartimento contém exatamente \(k\) pontos, então \(p_i ≈ k/n\) e, portanto, é razoável atribuir \[ h_i=\dfrac{k}{nf(x)}, \] o que implica \[ \tag{3.35} \mbox{AAMSE}(x)=\dfrac{f(x)^2}{k}+\dfrac{1}{12}\dfrac{k^2}{n^2}\left(\dfrac{f'(x)}{f(x)} \right)^2\cdot \] Agora, \(k\) deve ser escolhido não para minimizar o \(\mbox{AAMSE}(x)\) pontual, mas sim o global \[ \tag{3.36} \mbox{AAMISE}(k)=\dfrac{R(f)}{k}+\dfrac{1}{12}\dfrac{k^2}{n^2}\int \dfrac{f'(x)^2}{f(x)^2}\mbox{d}x\cdot \]

Agora \(k^∗ = O(n^{2/3})\) com \(\mbox{AAMISE}^∗ = O(n^{-2/3})\). Mas calcular a integral no termo de viés quando \(f =\phi\), a densidade normal padrão, dá \[ \tag{3.37} \int_\mathbb{R} \dfrac{f'(x)^2}{f(x)^2}\mbox{d}x = \int_{-\infty}^\infty \dfrac{\big(-x \, \phi(x) \big)^2}{\phi(x)^2}\mbox{d}x =\int_{-\infty}^\infty x^2 \mbox{d}x =\infty\cdot \]

Este resultado intrigante pode ser apenas consequência do uso da aproximação (3.35) combinada com o suporte infinito de \(\phi\), ou pode indicar que essas malhas sofrem vieses inesperadamente grandes.

Os resultados apresentados na Figura 3.15 demonstram claramente o desempenho insatisfatório das malhas percentuais quando comparadas com intervalos de largura fixa, com a diferença aumentando com o tamanho da amostra. As consequências, se houver, desse resultado negativo para os testes qui-quadrado são desconhecidas (ver Exercício 3.13).

Figura 3.15: MISE exato para malhas fixas (•••••) e percentis (ooooo) com \(k\) intervalos sobre (−1,1) para a densidade \(Beta(5,5)\) transformada com \(n = 100\) e 1000.


3.2.8.5 Uso de malhas adaptativas versus transformação


Em análises univariadas, as malhas adaptativas são mais úteis quando os dados são assimétricos ou quando estão agrupados com diferentes escalas dentro de cada grupo. Para dados com caudas pesadas, as malhas adaptativas podem proporcionar alguma melhoria visual, mas o critério MISE é em grande parte insensível a essas características. A melhor recomendação é usar a escada de transformações de Tukey (1977) ou os métodos de Box and Cox (1964) para reduzir a assimetria antes de aplicar uma malha fixa.

A escada de Tukey é um subconjunto da família de transformações de potência \(x^\lambda\), onde a transformação é definida como \(\log(x)\) quando \(\lambda = 0\): \[ \tag{3.38} \cdots, x^{-2},x^{-1},x^{-1/2},x^{-1/4},\log(x),x^{1/4},x^{1/2},x,x^2,x^3,\cdots \; \cdot \]

A família Box-Cox é semelhante, mas contínua em \(\lambda\) para \(x>0\): \[ \tag{3.39} x^{(\lambda)}=\left\{ \begin{array}{cc} \displaystyle \dfrac{x^\lambda-1}{\lambda}, & \lambda\neq 0 \\[0.8em] \log(x), & \lambda=0 \end{array}\right.\cdot \]

Transformações simples, como \(\sqrt{x}\) ou \(\log(x+c)\), geralmente são suficientes. No contexto multivariado, trabalhar com variáveis marginais fortemente assimétricas resulta em dados agrupados próximos às faces de um hipercubo circundante.

Mesmo transformações simples de variáveis marginais podem mover a nuvem de dados em direção ao centro do hipercubo, reduzindo significativamente um dos efeitos da maldição da dimensionalidade (ver Capítulo 7). No contexto univariado, Wand and Jones (1991) propuseram encontrar uma transformação ótima; dentro da família de transformações de Box-Cox, por exemplo, em relação ao AMISE nas coordenadas transformadas e, em seguida, transformar o estimador da densidade de volta para a escala original.

Em muitas aplicações, opta-se por trabalhar diretamente nas coordenadas transformadas. Economistas, por exemplo, geralmente trabalham com \(\log(\mbox{renda})\) em vez de renda para lidar com a forte assimetria presente nos dados de renda.


3.2.8.6 Observações


Histogramas adaptativos construídos de forma otimizada devem ser melhores do que histogramas não adaptativos de largura fixa. No entanto, histogramas adaptativos construídos de forma ad hoc não precisam ser melhores e, na verdade, podem ser muito piores.

Em particular, um algoritmo de histograma adaptativo pode ser inferior se não incluir larguras de intervalo fixas como um caso especial. Mesmo assim, na prática, um procedimento adaptativo pode ser vantajoso. Uma preocupação prática específica é a introdução de parâmetros de suavização adicionais que definem a adaptabilidade e que devem ser estimados a partir dos dados por algum método, como a validação cruzada (VC), que será apresentada posteriormente.

Dado o desempenho desses algoritmos no caso mais simples, com apenas um parâmetro de suavização em um cenário de largura de banda fixa, cautela parece ser necessária. Outros algoritmos de histograma adaptativo foram propostos por Van Ryzin (1973) e Wegman (1970).


3.3 Regras práticas de largura de intervalo amostrais


Embora ideias simples, como a regra de referência normal na equação (3.17), sejam úteis, a busca por procedimentos baseados em dados que minimizem aproximadamente o MISE e/ou o ISE é objeto de muita pesquisa. Os tópicos a seguir indicam a variedade desse trabalho.


3.3.1 Larguras de banda suavizadas em excesso


Em princípio, não há limite inferior para \(h\), pois a densidade desconhecida pode ser arbitrariamente irregular, embora para amostras finitas a largura de intervalo “ótima” possa ser surpreendentemente grande, como foi o caso para dados Lognormais.

Por outro lado, verifica-se que existem limites superiores úteis para a largura do intervalo, dependendo de algum conhecimento baseado em dados sobre a escala da densidade desconhecida (Terrell and Scott 1985).


3.3.1.1 Limites inferiores para o número de compartimentos


Reexaminando a expressão para o minimizador de AMISE na equação (3.15), \[ \tag{3.40} h^*=\left(\dfrac{6}{R(f')} \right)^{1/3} n^{-1/3}, \] conclui-se que qualquer limite inferior para \(R(f')\) leva a um limite superior para a largura do intervalo.

O conhecimento prévio mais simples da escala é que a densidade é zero fora de um intervalo \((a,b)\). Isso sugere o seguinte problema de otimização: \[ \tag{3.41} \min_f \int_{-\infty}^\infty f'(x)^2 \mbox{d}x, \qquad \mbox{suporte de } f = [-0.5,0.5]\cdot \]

Hodges and Lehmann (1956) resolveram esse problema de otimização em outro contexto e mostraram que a solução é \[ f_1(x)=\dfrac{3}{2}(1-4x^2)\pmb{I}_{-0.5,0.5}(x) = \dfrac{3}{2}(1-4x^2)_+\cdot \]

É instrutivo mostrar que \(f_1\) de fato minimiza \(R(f')\). Defina \(g(x) = f_1(x)+\epsilon(x)\), onde \(\epsilon(x)\) é uma função de perturbação que satisfaz as seguintes condições:

    1. \(\epsilon(x) = 0\) fora do intervalo \([-0.5,0.5]\);

    2. \(\epsilon(\pm 0.5) = 0\) e é contínua, pois caso contrário \(g(x)\) seria descontínua e \(R(g') =\infty\);

    3. \(\displaystyle \int \epsilon(x) \mbox{d}x = 0\), de modo que \(g\) integra para 1.

A solução deve ter suporte exatamente em \((-0.5, 0.5)\) e não em um subintervalo estrito. Caso contrário, reescalar a densidade para o intervalo \((-0.5, 0.5)\) reduziria a rugosidade, levando a uma contradição. Da mesma forma, a solução deve ser simétrica em torno de 0. Se \(g\) for assimétrico, considere a densidade simétrica \[ \dfrac{g(x) + g(-x)}{2}, \] cuja derivada é \(\big(g'(x)- g'(-x)\big)/2\).

Utilizando a desigualdade de Minkowski, que afirma que \[ \left( \int (f+g)^2\right)^{1/2} \leq \left(\int f^2 \right)^{1/2} + \left(\int g^2 \right)^{1/2}, \] a rugosidade da primeira derivada dessa densidade simétrica é igual a \[ \sqrt{\int \dfrac{1}{4} \Big(g'(x)-g'(-x) \Big)^2\mbox{d}x} \leq \dfrac{1}{2}\left( \sqrt{R(g')} + \sqrt{R(g'(-x)}\right)=\sqrt{R(g')}, \] o que contradiz a otimalidade de \(g\).

Para verificar a otimalidade de \(f_1\), calcule diretamente \[ \tag{3.42} \int g'(x)^2\mbox{d}x =\int f_1'(x)^2\mbox{d}x+\int \epsilon'(x)^2\mbox{d}x +2\int f_1'(x)\epsilon'(x)\mbox{d}x\cdot \] A última integral se anula quando \[ \int_{-1/2}^{1/2} f_1'(x)\epsilon'(x)\mbox{d}x =f_1'(x)\epsilon(x)\Bigg|_{x=-1/2}^{x=1/2}-\int_{-1/2}^{1/2} f''(x)\epsilon(x)\mbox{d}x=0, \] porque \(f_1''(x) = -12\), uma constante. Portanto, da equação (3.42), \(R(g') = R(f_1')+R(\epsilon')\), de modo que \(R(g')\geq R(f_1')\), o que prova o resultado. Observe que \(f_1\) é a única densidade quadrática que satisfaz todas as condições laterais.

Para aplicar esse resultado, observe que \(R(f') = 12\). Se \(f_1\) for mapeado linearmente para o intervalo \((a,b)\), então o limite inferior é \(R(f')\geq 12/(b-a)^3\). Portanto, \[ h^*=\left( \dfrac{6}{nR(f')}\right)^{1/3}\leq \left(\dfrac{6(b-a)^3}{n12} \right)^{1/3}=\dfrac{b-a}{\sqrt[3]{2n}}=h_{OS}, \] que é a chamada largura de banda excessivamente suavizada. Reorganizando, obtemos \[ \tag{3.43} \mbox{Número de intervalos } = \dfrac{b-a}{h^*}\geq \dfrac{b-a}{h_{OS}}=\sqrt[3]{2na}\cdot \]

Terrell and Scott (1985) mostraram que o intervalo da amostra pode ser usado se o intervalo \((a,b)\) for desconhecido ou mesmo se \(b-a = \infty\) e a cauda não for muito pesada.


Exemplo 3.2:

Para o conjunto de dados de queda de neve de Buffalo (\(n\) = 63), temos em snow a precipitação anual de neve em Buffalo, NY, entre 1910 e 1972, em polegadas.

snow = data.frame(snow = c(126.4,82.4,78.1,51.1,90.9,76.2,104.5,87.4,110.5,25.0,69.3,53.5,39.8,
         63.6,46.7,72.9,79.6,83.6,80.7,60.3,79.0,74.4,49.6,54.7,71.8,49.1,103.9,
         51.6,82.4,83.6,77.8,79.3,89.6,85.5,58.0,120.7,110.5,65.4,39.9,40.1,
         88.7,71.4,83.0,55.9,89.9,84.8,105.2,113.7,124.7,114.5,115.6,102.4,
         101.4,89.8,71.5,70.9,98.3,55.5,66.1,78.4,120.5,97.0,110.0))

Fonte: Carmichael (1976) e Parzen (1979).

Também dispomos dos dados LRL que contêm a contagem de intervalos de um experimento de física de partículas. Há 172 intervalos com centros variando de 285 a 1995 e largura de intervalo de 10 MeV. Um MeV é igual a um milhão de elétron-volts. O tamanho da amostra é 25.752.

A regra (3.43) sugere entre 5 e 37 intervalos, respectivamente. Para os dados de queda de neve, a faixa de amostragem é 126.4 - 25.0 \(\approx\) 100; portanto, sugere-se uma largura de banda suavizada de \(h_{OS} = 20\) polegadas de neve (ver Figura 3.16). Não há indícios do comportamento trimodal que parece estar presente no histograma com mais intervalos.

LRL = data.frame( LRL = c(5,11,17,21,15,17,23,25,30,22,36,29,33,43,54,55,59,44,58,
                          66,59,55,67,75,82,98,94,85,92,102,113,122,153,155,193,197,
                          207,258,305,332,318,378,457,540,592,646,773,787,783,695,
                          774,759,692,559,557,499,431,421,353,315,343,306,262,265,254,
                          225,246,225,196,150,118,114,99,121,106,112,122,120,126,126,
                          141,122,122,115,119,166,135,154,120,162,156,175,193,162,178,
                          201,214,230,216,229,214,197,170,181,183,144,114,120,132,109,
                          108,97,102,89,71,92,58,65,55,53,40,42,46,47,37,49,38,29,34,
                          42,45,42,40,59,42,35,41,35,48,41,47,49,37,40,33,33,37,29,26,
                          38,22,27,27,13,18,25,24,21,16,24,14,23,21,17,17,21,10,14,18,
                          16,21,6))

Fonte: Good and Gaskins (1980).

Os dados LRL foram registrados em 172 intervalos consecutivos, de modo que \[ h_{OS} = (787-5)/37 = 21.1\cdot \] Como os dados originais foram arredondados para intervalos de largura de 10 MeV, as únicas opções disponíveis para \(h_{OS}\) são 21 ou 22 MeV.

library(ggplot2)
library(patchwork)
p1 = ggplot(snow, aes(x = snow)) +
  geom_histogram(aes(y = ..density..), 
                 binwidth = 20,        
                 color = "white",       
                 fill = "#0073C2",      
                 alpha = 0.6) +         
  labs(title = "Buffalo, NY",
       x = "Precipitação anual de neve (em polegadas)",
       y = "Densidade") + # Renomeia os eixos e adiciona título
  theme_minimal() # Aplica um tema visual limpo

p2 = ggplot(LRL, aes(x = LRL)) +
  geom_histogram(aes(y = ..density..), 
                 binwidth = 22,        
                 color = "white",       
                 fill = "#0073C2",      
                 alpha = 0.6) +         
  labs(title = "LRL",
       x = "MeV",
       y = "Densidade") + 
  theme_minimal() 
p1 + p2

Figura 3.16: Histogramas suavizados em excesso dos dados de queda de neve (snow) em Buffalo e LRL.

Estritamente falando, combinar apenas quatro intervalos adjacentes \((h = 21)\) não é um alisamento excessivo, portanto, a escolha ligeiramente conservadora de \(h = 22\) é adotada (veja a Figura 3.16). Pelo menos três grupos aparecem neste histograma, o terceiro em torno de \(x = 200\).


3.3.1.2 Limites superiores para a largura de banda


O resultado anterior destaca-se pela sua simplicidade e facilidade de memorização. No entanto, o foco está na escolha do parâmetro de suavização e não no número de intervalos. A suavização excessiva também oferece uma solução nesse contexto.

Considere o problema de otimização (3.41), exceto que a restrição de intervalo fixo é substituída pela restrição de que a variância de \(f\) seja igual a \(\sigma^2\). Terrell (1990) mostrou que a solução para este problema é \[ f_2(x)=\dfrac{15}{16\sqrt{7} \sigma}\left(1-\dfrac{x^2}{7\sigma^2} \right)^2 \pmb{I}_{[-\sqrt{7}\sigma,\sqrt{7}\sigma]}(x)\cdot \]

Pode-se verificar que \(R(f'_2)=15\sqrt{7}/\big(343\sigma^3\big)\); portanto, \[ \tag{3.44} h^*\leq \left(\dfrac{6}{nR(f_2')} \right)^{1/3}=\left(\dfrac{686\, \sigma^3}{5\sqrt{7}n} \right)^{1/3}\approx 3.729 \, \sigma \, n^{-1/3}=h_{OS}\cdot \] Terrell (1990) considerou outras restrições de escala menos comuns. A versão baseada no intervalo interquartil (IQR) é particularmente robusta: \[ \tag{3.45} h^*\leq 2.603 (\mbox{IQR})n^{-1/3}=h_{OS}\cdot \]


Exemplo 3.3:

Para dados normais \(f=N(\mu,\sigma^2)\), do qual obtemos \[ \tag{3.46} h^*=2.5 \, \sigma \, n^{-1/3}< 3.729 \, \sigma \, n^{-1/3}=h_{OS}, \] na verdade, obtevemos o limite superior de \(h_{OS}\), mas \(h^*\) é apenas 93.6% de \(h_{OS}\).

Aparentemente, a densidade normal é muito suave e o uso da regra simples (3.17) é ainda mais justificado. No caso de dados Lognormais, para os quais \(\sigma^2=e(e-1)\), \[ \tag{3.47} h^*=1.44 \, n^{-1/3} < 3.729 \times 2.161 \times n^{-1/3}=8.059\, n^{-1/3}, \] que corresponde a 18% de \(h_{OS}\).

Por outro lado, para \(n = 25\), o melhor valor de \(h\approx 1.28\) da Figura 3.11, que corresponde a 46% de \(h_{OS} = 2.756\). Para dados de Cauchy, o limite baseado na variância é assintoticamente sem sentido, enquanto a regra de suavização excessiva baseada no intervalo interquartil (IQR) ainda se aplica.


3.3.2 Validação cruzada (VC) com e sem viés


Nesta seção, apresentamos a abordagem de validação cruzada (VC) para a calibração automática de histogramas baseada em dados. Os algoritmos de VC reutilizam os dados. O objetivo não é simplesmente produzir uma sequência consistente de larguras de intervalo — mesmo a regra de Sturges é consistente.

Em vez disso, o objetivo é produzir, de forma confiável, larguras de intervalo \(\widehat{h}_{CV}\) próximas de \(h^∗\) para um número finito de amostras ou, talvez de forma ainda mais otimista, larguras de intervalo que minimizem os erros ISE para amostras individuais.


3.3.2.1 Validação cruzada com viés


A única incógnita no AMISE é \(R(f')\), que pode ser estimada usando os dados disponíveis. Uma aproximação de \(f'(t_k)\) está disponível com base em uma diferença finita do histograma nos pontos médios dos intervalos \(B_k\) e \(B_{k+1}\), ou seja, \[ \widehat{f}'(t_k) = \frac{\nu_{k+1}/(nh)-\nu_k/(nh)}{h}\cdot \]

Uma estimativa potencial de \(R(f')\) é \[ \tag{3.38} \widehat{R}_1=\sum_k \Big(\widehat{f}'(t_k) \Big)^2 h=\dfrac{1}{n^2h^3}\sum_k \big(\nu_{k+1}-\nu_k\big)^2\cdot \] Pode-se demonstrar (ver Exercício 3.16) que \[ \tag{3.49} \mbox{E}(\widehat{R}_1)=R(f')+2/(nh^3)+O(h)\cdot \] Com suavização ótima, \(2/(nh^3)\) converge para \(R(f')/3\) pelo Teorema 3.3. Segue-se que \(\widehat{R}_1\) é um estimador enviesado de \(R(f')\), muito grande por um fator de um terço, de modo que \(\frac{3}{4}\widehat{R}_1\) é um estimador assintoticamente não enviesado de \(R(f')\).

Alternativamente, \[ \tag{3.50} \widehat{R}_h(f')=\dfrac{1}{n^2h^3} \sum_k \big(\nu_{k+1}-\nu_k\big)^2-\dfrac{2}{nh^3}\cdot \] Substituindo (3.50) na expressão AMISE (3.15), obtém-se uma estimativa de validação cruzada (CV) enviesada (BCV) do \(\mbox{MISE}(h)\): \[ \tag{3.51} \mbox{BCV}(h)=\dfrac{5}{6nh}+\dfrac{1}{12n^2h}\sum_k\big(\nu_{k+1}-\nu_k \big)^2, \] onde \(\nu_k\) é recalculado como \(h\) e a malha varia, para maior clareza, a origem da malha \(t_0\), permanece fixa.

O termo “bias” refere-se simplesmente ao fato de que o \(\mbox{AMISE}(h)\) é uma aproximação enviesada do verdadeiro \(\mbox{MISE}(h)\). A largura de banda do BCV \(\widehat{h}_{BCV}\), é definida como o minimizador de \(\mbox{BCV}(h)\), sujeito à restrição \(h\leq h_{OS}\): As propriedades teóricas dessas variáveis aleatórias foram estudadas por Scott and Terrell (1987), que mostraram que \(\widehat{h}_{BCV}\) converge para o \(h^∗\) com um erro relativo de \(O(n^{-1/6})\).

Essa taxa de convergência (relativa) é particularmente lenta. No entanto, o algoritmo parece útil na prática, especialmente para amostras grandes. Exemplos são apresentados na Seção 3.3.2.4.


3.3.2.2 Validação cruzada não viesada


Rudemo (1982) empreendeu a tarefa mais ambiciosa de tentar minimizar o erro \(L_2\) ISE real de um histograma para uma amostra individual.

Expandindo a expressão \[ \mbox{ISE}(h)=\int \Big(\widehat{f}(x)-f(x) \Big)^2 \mbox{d}x=R(\widehat{f})-2\int \widehat{f}(x)f(x)\mbox{d}x+R(f), \] Rudemo observou-se que o minimizador de \(\mbox{ISE}(h)\) não dependia da incógnita \(R(f)\) e que \(R(\widehat{f̂})\) podia ser facilmente calculado em forma fechada.

Além disso, ele observou que a segunda integral podia ser escrita como \(\mbox{E}\Big(\widehat{f}(X)\Big)\) onde a esperança é relativa ao ponto de avaliação e não à amostra aleatória \(x_1,\cdots,x_n\).

A metodologia de validação cruzada sugere a remoção de um ponto de dados e a utilização dos \(n-1\) pontos restantes para construir um estimador de \(\mbox{E}\Big(\widehat{f}(X)\Big)\). O \(n\)-ésimo ponto de dados é então avaliado na estimativa com o objetivo de determinar a qualidade do ajuste. Este passo é repetido \(n\) vezes, uma vez para cada ponto de dados, e os resultados são calculados em média.

Rudemo considerou o histograma \(\widehat{f}_{-i}(x)\), que é o histograma baseado nos \(n-1\) pontos da amostra, excluindo \(x_i\). Pode-se verificar que a variável aleatória observável \(\widehat{f}_{-i}(X_i)\) tem a mesma média que a variável aleatória não observável \(\mbox{E}\Big(\widehat{f}(X)\Big)\), embora baseada em uma amostra de \(n-1\) em vez de \(n\) pontos.

Ao excluir cada um dos \(n\) pontos de dados um de cada vez, Rudemo obteve uma estimativa estável de \(\mbox{E}\Big(\widehat{f}(X)\Big)\) calculando a média dos \(n\) casos. Com essa estimativa, ele propôs minimizar a função de validação cruzada por mínimos quadrados ou validação cruzada não viesada (\(\mbox{UCV}\)). \[ \tag{3.52} \mbox{UCV}(h)=R(f')-\dfrac{2}{n}\sum_{i=1}^n \widehat{f}_{-i}(x_i)\cdot \]

Para o histograma, os termos na \(\mbox{UCV}\) são facilmente avaliados. Por exemplo, \[ \widehat{f}_{-i}(x_i)=\dfrac{\nu_k-1}{(n-1)h}, \qquad \mbox{se} \qquad x_i\in B_k\cdot \] A fórmula computacional final para o \(\mbox{UCV}\) é \[ \tag{3.53} \mbox{UCV}(h)=\dfrac{2}{(n-1)h}-\dfrac{n+1}{n^2(n-1)h}\sum_k \nu_k^2\cdot \]

Veja o Exercício 3.17. A semelhança entre as fórmulas \(\mbox{UCV}\) e \(\mbox{BCV}\) é notável e justifica parcialmente a denominação BCV, visto que a BCV se baseia no AMISE em vez do MISE.

Métodos que se baseiam na fórmula AMISE são frequentemente denominados métodos plug-in (PI). Scott and Terrell (1987) também estudaram as propriedades distribucionais assintóticas das variáveis aleatórias UCV e \(\widehat{h}_{UCV}\).

Eles demonstraram que as larguras de banda UCV e BCV eram consistentes, mas que a convergência era lenta. Especificamente, \[ \sigma_{_{h_{CV}}}/h_{CV}=O(n^{-1/6})\cdot \]

Entretanto, Hall and Marron (1987c) e Hall and Marron (1987b) mostraram que esta é a mesma taxa de convergência do parâmetro de suavização que realmente minimiza o ISE. Hall and Marron (1987a) também estudaram a estimação de funcionais como \(R(f')\).


3.3.2.3 Problemas finais com BCV e UCV


Para uma amostra fixa, \(\mbox{BCV}(h)\) e \(\mbox{UCV}(h)\) exibem comportamentos não encontrados em \(\mbox{MISE}(h)\) para valores grandes e pequenos de \(h\), respectivamente. Por exemplo, quando \(h\to\infty\), todos os compartimentos, exceto um ou dois, estarão vazios.

Pode-se observar que \[ \tag{3.54} \lim_{t\to\infty} \mbox{BCV}(h)=\lim_{t\to\infty} \mbox{UCV}(h)=0\cdot \]

Claro que, \(\mbox{ISE}(h)\to R(f)\) quando \(h\to\infty\). Portanto, para \(\mbox{UCV}(h)\), o limite em (3.54) está correto, visto que o termo \(R(f)\) é omitido em sua definição. Contudo, para \(\mbox{BCV}(h)\), que é assintoticamente não negativo próximo a \(h^*\), o valor limite de 0 significa que o minimizador global está, na verdade, em \(h = \infty\).

Na prática, \(\widehat{h}_{BCV}\) é escolhido para ser um minimizador local com a restrição de ser menor que a largura de banda suavizada. Se não houver minimizador local dentro dessa região, a própria largura de banda suavizada é o minimizador com restrição.

Por outro lado, quando \(h\to 0\), as contagens de intervalos devem ser todas zero ou um, assumindo que os dados sejam contínuos (sem empates). Nesse caso, tanto \(\mbox{BCV}(h)\) quanto \(\mbox{UCV}(h)\) são aproximadamente iguais a \(1/(nh)\), a expressão correta para a variância integrada (IV). Observe que \(\mbox{VI}(h)\to\infty\) quando \(h\to 0\). No entanto, se os dados contiverem muitos empates, talvez devido a poucos dígitos significativos ou arredondamento, \(\mbox{UCV}(h)\) pode divergir para \(-\infty\).

Por exemplo, suponha que as \(n\) amostras consistam em \(n/m\) valores distintos, cada um replicado \(m\) vezes; então, \(\mbox{UCV}(h)\) é aproximadamente \((2-m)/(nh)\), que diverge para \(-\infty\) (em oposição a \(+\infty\)) se \(m > 2\) quando \(h\to 0\). Assim, o minimizador global de \(\mbox{UCV}(h)\) ocorre em \(h = 0\) nesse caso. Nessas situações, a divergência pode ocorrer apenas em uma pequena vizinhança em torno de \(h = 0\), e \(\widehat{h}_{UCV}\) é novamente definido como o mínimo local apropriado, se existir, mais próximo da largura de banda excessivamente suavizada; veja, no entanto, a Figura 3.19.


3.3.2.4 Aplicações


A Figura 3.17 exibe as curvas \(\mbox{BCV}\) e \(\mbox{UCV}\) para uma amostra \(N(0,1)\) de 1000 pontos. As escalas verticais são comparáveis, visto que a \(\mbox{UCV}\) é deslocada por uma constante fixa \(R(f)\). Assim, a função \(\mbox{UCV}\) acaba sendo muito mais ruidosa do que a função \(\mbox{BCV}\); entretanto, para outros dados mais complexos, o minimizador da função \(\mbox{BCV}\) pode ser bastante inadequado — tipicamente, ele é enviesado para larguras de intervalo maiores.

A função \(\mbox{UCV}\) pode ser ruidosa, mas seu minimizador é correto em média. Para esses dados, os dois minimizadores são \(h_{BCV} = 0.32\) e \(h_{UCV} = 0.36\), enquanto a largura do intervalo suavizada em excesso é \(h_{OS} = 0.37\), usando a regra da variância na equação (3.44). Escolher \(h_{UCV}\) dentre os muitos minimizadores locais fortes na função \(\mbox{UCV}\) é um problema prático.

Figura 3.17: Funções \(\mbox{BCV}\) e \(\mbox{UCV}\) para uma amostra \(N(0,1)\) de 1000 pontos.

Em conjunto, o sobressuavização, a validação cruzada bidimensional (\(\mbox{BCV}\)) e a validação cruzada unitária (\(\mbox{UCV}\)) constituem um conjunto poderoso de ferramentas para a escolha da largura do intervalo. Especificamente, as curvas de \(\mbox{UCV}\) e \(\mbox{BCV}\) devem ser mostradas graficamente em escalas comparáveis, marcando a localização do limite superior dado pela largura de banda sobressuavizada.

Recomenda-se um gráfico log-log, mas como a \(\mbox{UCV}\) é negativa, apenas \(\log(h)\) é mostrado no gráfico. O minimizador é localizado e a qualidade dessa validação cruzada é avaliada subjetivamente, examinando-se o quão bem articulado o minimizador parece. Observe falhas óbvias: ausência de minimizador local em \(\mbox{BCV}(h)\) ou a solução degenerada de \(\mbox{UCV}\) com \(h=0\).

Caso não haja consenso entre as soluções propostas, examine os histogramas correspondentes e escolha aquele que apresentar menor ruído local, especialmente próximo aos picos.

O segundo exemplo concentra-se na renda de 1983, amostrada de 5626 famílias alemãs (DIW 1983). Os dados foram transformados para escala logarítmica, como é comum entre economistas, a fim de reduzir a penalidade \(\mbox{MISE}\) devido à extrema assimetria. Quatro histogramas desses pontos de dados transformados são apresentados na Figura 3.18. Para mais detalhes, consulte Scott and Schmidt (1988), Scott and Schmidt (1989).

DIW = read.csv(file = "http://estatistica.c3sl.ufpr.br/~lucambio/Nonparam/DIW.csv")
par(mfrow=c(2,3),mar=c(2,2,1,1))
hist(log(DIW$x), breaks = 13, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "13 intervalos", col = "black", cex = 1.2) 
hist(log(DIW$x), breaks = 51, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "51 intervalos", col = "black", cex = 1.2)
hist(log(DIW$x), breaks = 68, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "68 intervalos", col = "black", cex = 1.2)
hist(log(DIW$x), breaks = 95, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "95 intervalos", col = "black", cex = 1.2)
hist(log(DIW$x), breaks = 200, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "200 intervalos", col = "black", cex = 1.2)
hist(log(DIW$x), breaks = 300, col = "lightblue", main = "", freq = FALSE, ylim = c(0.0,0.8))
box();grid();text(x = 7, y = 0.5, labels = "300 intervalos", col = "black", cex = 1.2)

Figura 3.18: Seis histogramas da amostra de renda alemã de 1983, composta por 5625 famílias.

Qual histograma é o melhor? A regra de Sturges sugere apenas 13 intervalos, o que é claramente um suavizamento excessivo. A regra de suavização excessiva (largura máxima do intervalo) resulta em \(h = 0.160\), correspondendo a 68 intervalos na faixa da amostra. As funções \(\mbox{BCV}\) e \(\mbox{UCV}\) são apresentadas na Figura 3.19.

O minimizador do critério \(\mbox{BCV}\) ocorre em \(h = 0.115\), sugerindo 95 intervalos. A função \(\mbox{UCV}\) parece ter seu mínimo global em \(h = 0\), mas não apresenta um mínimo local significativo na região de interesse. De fato, a análise dos dados originais revelou apenas 3322 rendimentos únicos, portanto, o problema de dados duplicados está claramente presente. Mostramos a continuação a dimensão do arquivo de dados DIW e como obter extrair o número de elementos únicos, para isto utilizamos a função unique(), a qual retorna um vetor, um arquivo de dados ou um arranjo com os elementos/linhas duplicados removidos.

dim(DIW)
## [1] 5625    1
dim(unique(DIW))
## [1] 3322    1


Foi determinado empiricamente que esse efeito poderia ser eliminado adicionando-se um pequeno ruído uniforme aos valores logarítmicos, como mostrado na Figura 3.19. No entanto, ainda existe uma ampla gama de mínimos possíveis. Mesmo a pequena largura de intervalo \(h = 0.055\) parece viável, o que corresponde a 200 intervalos. Quando o ruído \(U(-0.005, 0.005)\) foi adicionado, a curva \(\mbox{UCV}\) tornou-se horizontal, dividindo as duas curvas mostradas na figura.

Figura 3.19: Funções \(\mbox{BCV}\) e \(\mbox{UCV}\) da amostra de renda alemã de 1983, composta por 5625 domicílios. Uma segunda curva \(\mbox{UCV}\) é mostrada para os dados logarítmicos suavizados.

Para grandes quantidades de ruído aditivo, a curva \(\mbox{UCV}\) foi visualmente afetada para larguras de banda próximas de \(h_{OS}\). Portanto, há um certo grau de subjetividade no processo de remoção de dados duplicados.

A lição é que todos os três algoritmos devem ser examinados simultaneamente, mesmo com amostras muito grandes. Este tópico será abordado com mais detalhes no Capítulo 6.


3.4 Teoria \(L_2\) para histogramas multivariados


A derivação do \(\mbox{MISE}\) para o histograma multivariado é apenas ligeiramente mais complicada do que no caso univariado.

Dado uma amostra de \(f(\pmb{x})\), onde \(\pmb{x}\in\mathbb{R}^d\), o histograma é determinado por uma partição do espaço. Considere uma partição regular por hiper-retângulos de tamanho \(h_1\times h_2\times \cdots\times h_d\). Escolher hipercubos como compartimentos seria suficiente se os dados estivessem devidamente escalonados, mas, em geral, esse não será o caso.

Melhorias adicionais podem ser obtidas considerando compartimentos não regulares ou rotacionados, ver Scott (1988) e Hüsemann and Terrell (1991).

Considere um recipiente hiper-retangular genérico rotulado como \(B_k\) contendo \(\nu_k\) pontos. Como de costume, \(\sum_k \nu_k=n\). Então \[ \widehat{f}(\pmb{x})=\dfrac{\nu_k}{n\, h_1h_2\cdots h_d} \qquad \mbox{para} \qquad \pmb{x}\in B_k\cdot \] A variância do histograma é constante em cada intervalo e é dada por \[ \tag{3.55} \mbox{Var}\big(\widehat{f}(\pmb{x})\big)=\dfrac{n\, p_k(1-p_k)}{\big(n\, h_1h_2\cdots h_d\big)^2} \qquad \mbox{para} \qquad \pmb{x}\in B_k\cdot \]

Integrando a variância ao longo do intervalo, multiplicamos (3.55) pelo volume do hiper-retângulo, \(h_1 h_2\cdots h_d\). Somando as contribuições da variância de todos os intervalos, obtemos a variância integrada como \[ \tag{3.56} \mbox{IV}=\dfrac{1}{n\, h_1h_2\cdots h_d}-\dfrac{R(f)}{n}+o(n^{-1}), \] onde \(R(f)=\displaystyle \int_{\mathbb{R}^d} f(\pmb{x})^2\mbox{d}\pmb{x}\).

O primeiro termo em (3.56) é exato, uma vez que \(\displaystyle \sum_k p_k = 1\); o resto é obtido observando que \[ p_k^2 \approx \big(h_1 h_2 \cdots h_d \, f(\epsilon_k)\big)^2, \] onde \(\epsilon_k\in B_k\) e usando aproximações de integração Riemanniana multivariada padrão.

O seguinte esboço para o cálculo do viés pode ser tornado rigorosamente. Considere o intervalo \(B_0\) centrado na origem \(\pmb{x} = \pmb{0}_d\). Agora \[ f(\pmb{x})=f(\pmb{0})+\sum_{i=1}^d x_i f_i(0)+\dfrac{1}{2}\sum_{i=1}^d\sum_{j=1}^d x_i x_j f_{ij}(\pmb{0})+O(h^3), \] onde \(f_i(\pmb{x})=\partial f(\pmb{x})/\partial x_i\) e \(f_{ij}(\pmb{x})=\partial^2 f(\pmb{x})/\partial x_i\partial x_j\). Portanto, \[ \tag{3.57} p_0=\int_{-h_d/2}^{h_d/2} \cdots \int_{-h_d/2}^{h_d/2} f(\pmb{x})\mbox{d}\pmb{x}= h_1 h_2\cdots h_d f(\pmb{0})+O(h^{d+2}), \] onde \(h=\min_i (h_i)\). Calculando o \(\mbox{Viés}\big(\widehat{f}(\pmb{x})\big)\) obtemos \[ \tag{3.58} \mbox{E}\big(\widehat{f}(\pmb{x})\big)-f(\pmb{x})=\dfrac{p_0}{n\, h_1 h_2\cdots h_d}-f(\pmb{x})=-\sum_{i=1}^d x_i f_i(\pmb{0})+O(h^2)\cdot \] Elevando ao quadrado e integrando sobre \(B_0\), o viés quadrático integrado para o intervalo é \[ \tag{3.59} h_1 h_2\cdots h_d \left( \sum_{i=1}^d \dfrac{1}{12}h_i^2 f_i(\pmb{0})^2 +O(h^4)\right)\cdot \] Uma expressão semelhante vale para todos os outros intervalos, com a origem \(\pmb{x} = \pmb{0}\) substituída pelos respectivos centros dos intervalos multivariados.

A soma de todos os intervalos resulta em \[ \tag{3.60} \mbox{AISB}(\pmb{h})=\dfrac{1}{12}\sum_{i=1}^d h_i^2 R(f_i)\cdot \] Essas aproximações estão reunidas em um teorema.


Teorema 3.5:

Para uma função de densidade suficientemente suave \(f(x)\), o \(\mbox{MISE}\) multivariado é assintoticamente \[ \tag{3.61} \mbox{AMISE}(\pmb{h})=\mbox{AIV}+\mbox{AISB}=\dfrac{1}{n\, h_1 h_2\cdots h_d}+\dfrac{1}{12}\sum_{i=1}^d h_i^2 R(f_i)\cdot \] As larguras de intervalo assintoticamente ótimas \(h^∗_k\) e o \(\mbox{AMISE}^∗\) resultante são \[ \tag{3.62} h_k^*=R(f_k)^{-1/2}\left(6\displaystyle \prod_{i=1}^d R(f_i)^{1/2} \right)^{1/(2+d)}n^{-1/(2+d)}, \] e \[ \tag{3.63} \mbox{AMISE}^*=\dfrac{d+2}{2}6^{-d/(2+d)}\left(\prod_{i=1}^d R(f_i) \right)^{1/(2+d)} n^{-2/(2+d)}\cdot \]


Demonstração ver Scott (1988) e Hüsemann and Terrell (1991).


Exemplo 3.?:

Suponha que \(X\sim N(\mu,\Sigma)\), \(\Sigma = \mbox{diag}(\sigma_1^2,\sigma_2^2,\cdots,\sigma_d^2)\). Então, considerando \(c_d=\sigma_1\sigma_2\cdots \sigma_d\), temos \[ \tag{3.64} R(f_i)=\big(2^{d+1}\pi^{d/2}\sigma_i^2 c_d \big)^{-1} \] e \[ h_k^*=2\, 3^{1/(2+d)} \pi^{d/(4+2d)}\sigma_k n^{-1/(2+d)}, \] \[ \tag{3.65} \mbox{AMISE}^*=(2+d) 2^{-(1+d)} 3^{-d/(2+d)}\pi^{d^2/(4+2d)}c_d^{-1}n^{-2/(2+d)}\cdot \]

Observe que a constante na largura de banda aumenta rapidamente de 3.4908 em uma dimensão até o valor limite de \(2\sqrt{\pi} = 3.5449\) quando \(d\to\infty\). Portanto, uma fórmula muito útil para memorizar é \[ \tag{3.66} \mbox{Regra de referência normal: } \; h_k^*\approx 3.5 \, \sigma_k n^{-1/(2+d)}\cdot \]


3.4.1 Maldição da dimensionalidade


Bellman (1961) cunhou a expressão “maldição da dimensionalidade” para descrever o crescimento exponencial na otimização combinatória à medida que a dimensão aumenta. Aqui, é o número de compartimentos que cresce exponencialmente com o aumento da dimensão. É importante observar como os histogramas são afetados por esse fenômeno.

Uma densidade relativamente simples de estimar é a normal multivariada com \(\Sigma=\pmb{I}_d\). Dando continuidade aos resultados do exemplo anterior, temos a Tabela 3.6. Claramente, a taxa de diminuição do \(\mbox{MISE}\) em relação ao tamanho da amostra degrada-se rapidamente à medida que a dimensão aumenta, em comparação com a taxa paramétrica ideal \(O(n^{-1})\). Uma grande vantagem da modelagem paramétrica é que a taxa de diminuição do \(\mbox{MISE}\) é independente da dimensão.

Tabela 3.6: Exemplo de larguras de intervalo e erros assintoticamente ótimos para \(f=N(\pmb{0}_d,\pmb{I}_d)\).

Entretanto, diversos analistas de dados têm utilizado histogramas na classificação de dados quadrivariados de sensoriamento remoto com resultados satisfatórios (Wharton 1983). Uma possível observação otimista é que as constantes do \(\mbox{AMISE}\) na Tabela 3.6 também diminuem à medida que a dimensão aumenta. Infelizmente, o \(\mbox{MISE}\) não é uma grandeza adimensional e, portanto, esses coeficientes não são diretamente comparáveis.

Epanechnikov (1969) descreveu um procedimento para comparar erros e desempenho de histogramas em diferentes dimensões. Ele considerou uma possível reescala adimensional do \(\mbox{MISE}\): \[ \tag{3.67} \epsilon_d=\dfrac{\mbox{MISE}}{R(f)} \approx \dfrac{2+d}{2}3^{-\frac{d}{2+d}}\pi^{\frac{d}{2+d}}n^{-\frac{2}{2+d}}, \qquad \mbox{quando} \qquad f=N(\pmb{0}_d,\pmb{I}_d)\cdot \]

Novamente, para o caso normal, pode-se calcular uma tabela de tamanhos de amostra equivalentes (ver Tabela 3.7). Esta tabela ilustra graficamente a maldição da dimensionalidade e a intuição de que a estimativa de densidade em mais de duas ou três dimensões não funcionará. Esta conclusão é, no entanto, muito pessimista, e outras evidências serão apresentadas na Seção 7.2.

Tabela 3.7: Tamanhos de amostra equivalentes entre dimensões para a densidade normal multivariada, com base no critério de Epanechnikov.

Outra abordagem para entender esse problema é contar os compartimentos em uma região de interesse. Para este mesmo exemplo normal, considere a região de interesse como uma esfera de raio \(r_d\), contendo 99% da massa de probabilidade.

O raio \(r_d\) é a solução para \[ \tag{3.68} P\left( \sum_{i=1}^d Z_i^2 \leq r_d^2\right) = 0.99, \qquad \mbox{logo} \qquad r_d=\sqrt{\chi^2_{0.99}}, \]

Ou seja, a raiz quadrada do quantil qui-quadrado apropriado. O volume da esfera é dado pela equação (1.3). O número de compartimentos do hipercubo sobre os dados nesta esfera é aproximadamente igual ao volume da esfera dividido pelo volume do compartimento \(h^d\), mais os compartimentos que cobrem 1% dos pontos fora da esfera.

A Tabela 3.8 ilustra esses cálculos para uma amostra de tamanho 1000. Em cinco dimensões, o histograma suavizado otimamente tem aproximadamente 1250 compartimentos na região de interesse.

Tabela 3.8: Número aproximado de intervalos na região de interesse para um histograma multivariado de 1000 pontos normais.

Como existem apenas 1000 pontos, é evidente que a maioria desses compartimentos estará vazia e que o histograma será bastante irregular. Scott and Thompson (1983) denominaram isso de “fenômeno do espaço vazio”. À medida que a dimensão aumenta, o histograma fornece estimativas razoáveis apenas perto da moda e longe das caudas.


3.4.2 Um caso especial: \(d = 2\) com correlação não nula


Os efeitos da correlação, se houver, foram ignorados até este ponto. Considere o caso normal bivariado simples, \(f(x_1,x_2)=N(\mu_1,\mu_2,\sigma_1^2,\sigma_2^2,\rho)\), para a qual \[ R(f_1)=\dfrac{1}{8\pi (1-\rho^2)^{3/2}\sigma_1^3\sigma_2} \] e \[ R(f_2)=\dfrac{1}{8\pi(1-\rho^2)^{3/2}\sigma_1\sigma_2^3}\cdot \] Do Teorema 3.5, \[ h^*=3.504 \, \sigma_1 (1-\rho^2)^{3/8} n^{ -1/4} \] e \[ \tag{3.69} \mbox{AMISE}^*=\dfrac{0.122}{\sigma_1\sigma_2} (1-\rho^2)^{-3/4}n^{-1/2}\cdot \]

O efeito da correlação \(\rho\) é introduzir potências da quantidade \((1-\rho^2)\) nas equações. Assim, se os dados não forem independentes, mas estiverem agrupados ao longo de uma linha, larguras de intervalo menores serão necessárias para “rastrear” essa característica. Se a densidade for degenerada, ou seja, \(\rho = \pm 1)\), o \(\mbox{MISE}\) explode. Esse resultado também indica que, se os dados caírem em qualquer variedade de dimensão inferior, um histograma nunca será consistente!

Essa inconsistência é um segundo aspecto, talvez mais importante, da “maldição da dimensionalidade” aplicada à estatística. Portanto, uma parte significativa do esforço em uma boa análise de dados deve ser dedicada a verificar a classificação (local) dos dados e identificar estruturas não lineares em dimensões inferiores (ver Capítulo 7).


3.4.3 Malhas bivariadas regulares ótimas


Considere o mosaico do plano com polígonos regulares como compartimentos para um histograma bivariado. Além dos quadrados, existem apenas duas outras malhas poligonais regulares: triângulos equiláteros e hexágonos (ver Figura 3.20). O mosaico não regular do plano é, obviamente, viável, mas é de interesse secundário.

Figura 3.20: As três possíveis malhas poligonais regulares para histogramas bivariados.

Scott (1988) comparou o \(\mbox{AMISE}\) bivariado para essas três formas de mosaico regulares. Se os compartimentos individuais forem parametrizados de forma que cada um tenha a mesma área \(h^2\), então Scott mostrou que \[ \tag{3.70} \mbox{AMISE}(h)=\dfrac{1}{nh^2}+ch^2\Big(R(f_x)+R(f_y)\Big), \] onde \(f_x=\partial f(x,y)/\partial x\) e \(f_y=\partial f(x,y)/\partial y\), e \[ \tag{3.71} c=\left(\dfrac{1}{12},\dfrac{1}{6\sqrt{3}},\dfrac{5}{36\sqrt{3}} \right)=\left(\dfrac{1}{12},\dfrac{1}{10.39},\dfrac{1}{12.47} \right), \] para os ladrilhos quadrados, triangulares e hexagonais, respectivamente (ver Exercício 3.25).

Portanto, as malhas hexagonais são de fato as melhores, mas apenas marginalmente. As malhas triangulares são bastante inferiores. Scott também considerou malhas compostas por triângulos retângulos. Ele observou que estas eram inferiores aos triângulos equiláteros e que o \(\mbox{ISB}\) às vezes incluía termos de produto vetorial, como \[ \int\int f_x(x,y)f_y(x,y) \mbox{d}x\mbox{d}y\cdot \]

Carr et al. (1987) sugeriram o uso de compartimentos hexagonais por um motivo diferente. Eles estavam avaliando a plotagem de glifos bivariados com base em contagens em compartimentos bivariados. O uso de compartimentos quadrados resultou em um alinhamento visualmente perturbador dos glifos nas direções vertical e horizontal. Essa perturbação foi praticamente eliminada quando os glifos foram posicionados de acordo com uma malha hexagonal. A eliminação é facilmente perceptível nos exemplos da Figura 3.21.

Figura 3.21: Glifos hexagonais para os dados de duração defasada do Old Faithful e para o conjunto de dados lipídicos de 320 homens com doença cardíaca.

Observe que o tamanho de cada glifo hexagonal é desenhado de forma que sua área seja proporcional ao número de observações naquele compartimento. Os autores observam que o método é especialmente útil para conjuntos de dados muito grandes. Observe que, embora o critério \(\mbox{MISE}\) não diferencie fortemente entre compartimentos quadrados ou hexagonais, fatores subjetivos podem influenciar fortemente a equação, para um lado ou para o outro.


3.5 Modas o picos em um histograma


Muitos dos exemplos apresentados até aqui mostraram dados multimodais. Good and Gaskins (1980) discutiram métodos avançados para encontrar modas e picos em dados. Em \(1-\mbox{D}\), uma moda é um conjunto (uma coleção de pontos) onde \(f'(x) = 0\) e \(f''(x) < 0\). Um pico é um conjunto (uma coleção de intervalos disjuntos) onde \(f''(x) < 0\). Assim, a busca por picos é mais geral do que a estimativa da moda. A Figura 3.22 mostra um histograma que contém uma moda e um pico. Um pico não contém necessariamente uma moda, embora uma moda esteja sempre localizada em um pico.

Espessura em milésimos de polegada de uma amostra de 90 moedas de um centavo de dólar dos EUA, cunhadas entre 1945 e 1989. Duas moedas foram medidas para cada ano. 1 milésimo de polegada = 0.001 polegada.

penny = read.csv(file = "http://estatistica.c3sl.ufpr.br/~lucambio/Nonparam/penny.csv", 
                 header = TRUE, sep = ",")
par(mfrow=c(1,2),mar=c(2,2,1,1))
hist(penny$mm, breaks = 10, col = "lightblue", main = "", freq = FALSE)
box();grid()
hist(LRL$LRL, breaks = 10, col = "lightblue", main = "", freq = TRUE)
box();grid()

m1_pois <- glm(LRL ~ 1, data = LRL, family = gaussian)
topmodels::rootogram(m1_pois, plot = "ggplot2",style = "standing")

Figura 3.22: Histograma da espessura de uma moeda de um centavo de dólar americano com uma moda e um pico. No quadro à direita, o histograma do conjunto de dados LRL com \(h = 10\) MeV é mostrado em escala de raiz quadrada; os 13 picos encontradas por Good e Gaskins são indicadas pelos segmentos de linha acima do histograma.

Em situações práticas de busca por picos, a densidade é frequentemente considerada uma mistura de várias densidades componentes (veja Izenman and Sommer 1988, por exemplo). O problema da mistura é ainda mais geral do que a busca por picos, já que uma densidade de mistura normal como \[ \tag{3.72} f(x) = \sum_{i=1}^q \omega_i \phi(x \, | \, \mu_i,\sigma_i^2) \qquad \mbox{onde} \qquad \sum_{i=1}^q \omega_i=1, \] não precisam exibir modos ou picos (ver Exercício 3.26).

No entanto, a estimativa dos parâmetros \(\{q, \omega_i,\mu_i,\sigma_i^2\}\) é frequentemente mal condicionada, especialmente se \(q\) for maior que o número verdadeiro de densidades Redner and Walker (1984).

Com dados de histograma, é natural examinar gráficos de primeiras e segundas diferenças (padronizadas) em busca de evidências de modas e picos, respectivamente: \[ \dfrac{\nu_{k+1}-\nu_k}{\sqrt{\nu_{k+1}+\nu_k}} \quad \mbox{e} \quad \dfrac{\nu_{k+1}-2\nu_k+\nu_{k-1}}{\sqrt{\nu_{k+1}+4\nu_k+\nu_{k-1}}}, \] visto que as contagens nos intervalos podem ser aproximadamente modeladas como variáveis aleatórias de Poisson independentes. Por exemplo, \(\mbox{Var}(\nu_{k+1}-\nu_k)\approx \nu_{k+1}+\nu_k\).

Figura 3.23: Diferenças de segunda ordem da raiz quadrada do conjunto de dados do histograma LRL com largura de intervalo de 30 MeV. Os 13 picos encontrados por Good e Gaskins são indicados como anteriormente. A linha tracejada indica o nível de corte aproximado de 5% para que um pico seja considerado significativo.

Alternativamente, podemos nos concentrar nas diferenças finitas do rootgrama, usando o fato de que a raiz quadrada é a transformação estabilizadora de variância para variáveis aleatórias de Poisson, com variância de 1/4: \[ \tag{3.73} \sqrt{\nu_{k+1}}-\sqrt{\nu_k} \qquad \mbox{e} \qquad \sqrt{\nu_{k+1}}-2\sqrt{\nu_k}+\sqrt{\nu_{k-1}}; \] que possuem variâncias aproximadas de \(1/2 = \frac{1}{4}+\frac{1}{4}\) e \(3/2 = \frac{1}{4} +(-2)^2\frac{1}{4}+\frac{1}{4}\), respectivamente, assumindo independência.

Good e Gaskins encontraram 13 picos em seu conjunto de dados LRL, veja a Figura 3.22. Um gráfico das segundas diferenças das contagens dos intervalos de raiz com \(h = 30\) MeV é mostrado na Figura 3.23. Sob a hipótese nula de que não há pico, ou seja, a segunda diferença é não negativa, o nível de significância de 5% para o teste unilateral será de −1.645 = −2.015$, o que é indicado pela linha tracejada na Figura 3.23.

Lembre-se de que esses dados foram coletados em intervalos de largura de 10 MeV. Para intervalos menores, a figura apresentava muito ruído. Ambos os critérios de validação cruzada indicam \(h = 10\) MeV como a melhor largura de intervalo. A maioria das protuberâncias encontradas por Good e Gaskins pode ser vista no gráfico. No entanto, algumas protuberâncias não são facilmente identificáveis. Uma grande porção negativa da curva próxima a \(x = 1145\) MeV (indicada pelo triângulo sólido) parece ter passado despercebida.


3.5.1 Propriedades do histograma “Modas”


Suponha, sem perda de generalidade, que a densidade verdadeira \(f\) tenha uma moda em 0. Considere um histograma com espaçamento uniforme, com o intervalo \(B_0\) centrado em \(x = 0\). A contagem de intervalos \(\nu_0\sim B(n,p_0 )\approx P(\lambda_0 )\), que é a densidade de Poisson com \(\lambda_0 = n p_0\). Assintoticamente, as contagens de intervalos adjacentes, \((\nu_{-k},\cdots,\nu_0,\cdots,\nu_k)\), são independentes e normalmente distribuídas com \(\nu_i ≈ N(\lambda_i,\lambda_i)\).

Considere a seguinte questão: Qual é a probabilidade de o histograma ter uma moda amostral no intervalo \(B_0\)? Condicionando à contagem de intervalos observada em \(B_0\), obtemos \[ \tag{3.74} \begin{array}{rcl} \displaystyle P\Big( \nu_0=\underset{{|j|\leq k}}{\mbox{arg max} \; \nu_j} \Big) & = & \displaystyle \sum_{x_{_0}} P\Big( \nu_0=\underset{{|j|\leq k}}{\mbox{arg max} \; \nu_j} \, \Big| \, \nu_0=x_0\Big)f_{\nu_{_{0}}}(x_0) \\[0.8em] & = & \displaystyle P\Big( \nu_j<x_0; \, |j|\leq k, \, j\neq 0 \Big)f_{\nu_{_{0}}}(x_0) \\[0.8em] & \approx & \displaystyle \int_{x_{_{0}}} \prod_{j=-k \\ j\neq 0}^k \Phi\left(\dfrac{x_0-\lambda_j}{\sqrt{\lambda_j}}\right) \phi\left(\dfrac{x_0-\lambda_0}{\sqrt{\lambda_0}} \right)\dfrac{1}{\sqrt{\lambda_0}}\mbox{d}x_0, \end{array} \] usando a aproximação normal para \(P(\nu_j < x_0)\), uma aproximação normal para a densidade de Poisson discreta \(f_{\nu_{0}}(x_0)\) e substituindo a soma por uma integral. Agora \(\lambda_j = n p_j\).

Observando que \(f'_0 ≡ f'(0) = 0\), temos como aproximação para \(p_j\): \[ \tag{3.75} p_j = \int_{(j-1/2)h}^{(j+1/2)h} \left(f_0+\dfrac{x^2}{2}f''_0+\cdots \right)\mbox{d}x = hf_0+\dfrac{h^3}{6}\left(3j^2+\dfrac{1}{4} \right)f''_0+\cdots \, \cdot \]

Realizando a mudança de variáveis \(y=(x_0-\lambda_0)/\sqrt{\lambda_0}\), segue-se que (3.74) é igual a \[ \tag{3.76} \begin{array}{rcl} \displaystyle P\Big( \nu_0=\underset{{|j|\leq k}}{\mbox{arg max} \; \nu_j} \Big) & \approx & \displaystyle \int_{y} \prod_{j=-k \\ j\neq 0}^k \Phi\left(\dfrac{\lambda_0-\lambda_j+y\sqrt{\lambda_0}}{\sqrt{\lambda_j}}\right) \phi\left(y\right)\mbox{d}y\\[0.8em] & \approx & \displaystyle \int_{y} \prod_{j=-k \\ j\neq 0}^k \Phi\left(y-\dfrac{j^2h^{5/2}\sqrt{n}}{2}\dfrac{f''(0)}{\sqrt{f(0)}}\right) \phi\left(y\right)\mbox{d}y\cdot \end{array} \]

No caso de um histograma ótimo, \(h=c\, n^{-1/3}\), do qual, \[ \lim_{n\to\infty} \Big(h^{5/2}\sqrt{n} \Big)=O\big(n^{-1/3}\big) =0, \] portanto, \[ \tag{3.77} \lim P\Big( \nu_0=\underset{{|j|\leq k}}{\mbox{arg max} \; \nu_j}\Big) = \displaystyle \int_y \Phi(y)^{2k}\phi(y)\mbox{d}y=\dfrac{1}{2k+1}, \] da equação (3.76).

A probabilidade \(1/(2k+1)\) está longe de um valor mais desejável próximo de 1. A interpretação correta desse resultado é que o \(\mbox{MISE}\) ótimo alisado resulta em intervalos com larguras muito estreitas para estimar as modas. De fato, em uma vizinhança da moda, a densidade parece plana e os intervalos têm essencialmente a mesma altura esperada em primeira ordem, de modo que cada um dos \(2k+1\) intervalos tem a mesma probabilidade de ser a moda amostral.

Em uma simulação de dados normais suavizados otimamente com \(n = 256000\), o número médio de modas amostrais no histograma foi 20! É verdade que a maioria dessas “modas” eram apenas pequenas aberrações, mas o resultado é bastante inesperado, veja a Figura 3.24.

Figura 3.24: Diagrama de raízes de dois histogramas de um milhão de pontos normais.

Em seguida, suponha que a origem não seja uma moda; então \(f'(0) = 0\). Uma análise semelhante mostra que a probabilidade de o intervalo \(B_0\) ser uma moda, o que não é!, converge para uma probabilidade fixa não nula quando \(n\to\infty\). A probabilidade é menor quanto maior for a magnitude de \(f'(0)\). Se intervalos mais amplos forem usados, as probabilidades podem convergir para os valores desejados de 1 e 0, respectivamente.

Um caso especial interessante para o limite na equação (3.77) decorre do uso de uma largura de intervalo muito maior \(h=c\, n^{-1/5}\). Essa escolha será examinada no Capítulo 4, que trata do estimador de densidade de polígonos de frequência intimamente relacionado.


3.5.2 Ruído em histogramas ótimos


Considere novamente os histogramas de um milhão de pontos normais mostrados na Figura 3.7. Ao mostrar um rootgrama, que é o histograma com variância estabilizada mostrado em escala de raiz quadrada, para os casos \(h = h^∗\) e \(h = h^∗/2\), fica evidente a presença de muitos modos falsos locais, veja a Figura 3.24. A moda dominante em zero, no entanto, não se perde em meio a todos os modos ruidosos.

Em resumo, histogramas que são “ótimos” em relação ao \(\mbox{MISE}\) podem não ser “ótimos” para outros fins. Essa constatação específica sobre as modas não deve ser surpreendente, dada a presença de ruído assintótico fixo em \(R(\widehat{f}')\) encontrado na derivação de validação cruzada enviesada. Mas essa característica geral dos procedimentos não paramétricos é bastante diferente daquela observada no contexto paramétrico, em que uma estimativa de densidade de máxima verossimilhança será ótima para praticamente qualquer aplicação específica.

No contexto não paramétrico, a calibração “ótima” dependerá da finalidade. Para amostras pequenas, o ruído no histograma raramente é atribuído à fonte correta e, portanto, não é bem compreendido. De qualquer forma, histogramas suavizados otimamente não fornecem uma ferramenta poderosa para a detecção de picos, veja novamente o exemplo LRL. Na verdade, as segundas diferenças de um histograma suavizado otimamente divergem da verdadeira segunda derivada. Esse resultado será mais fácil de visualizar posteriormente com estimadores que sejam diferenciáveis.


3.5.3 Larguras de banda ótimas do histograma para modas


Nesta seção, buscamos uma largura de banda que forneça melhores estimativas da derivada da densidade desconhecida. Usando a notação de intervalos como na Figura 3.2, uma estimativa natural da primeira derivada \(f'(x)\) é a primeira diferença finita \[ \widehat{f}'(x)=\dfrac{\widehat{f}_0(x)-\widehat{f}_{-1}(x)}{h}=\dfrac{\nu_0-\nu_{-1}}{nh^2}, \qquad -h/2<x<h/2\cdot \]

A série de Taylor de \(f(x)\) em torno de \(x = 0\) leva às aproximações \[ hf(0)\pm h^2\dfrac{f'(0)}{2}+h^3\dfrac{f''(0)}{6}\pm h^4\dfrac{f'''(0)}{24}+O(h^5) \] para \(p_0\) e \(p_{-1}\). Dado que \(\mbox{E}(\nu_i)=np_i\), \[ \mbox{E}\big(\widehat{f}'(x) \big)=\dfrac{np_0-np_{-1}}{nh^2}=f'(0)+\dfrac{1}{12}h^2 f'''(0)+\cdots, \qquad -h/2<x<h/2\cdot \]

Agora \(f'(x)=f'(0)+xf''(0)+x^2f'''(0)/2+\cdots\), então o viés é \(-xf''(0)+O(h^2)\). Então, o viés quadrático integrado sobre o intervalo é \(f''^2(0) h^2/12\times h\). Substituindo 0 por \(kh\) no \(k\)-ésimo intervalo e somando, obtemos \[ \mbox{AISB}=h^2\dfrac{R(f''')}{12}\cdot \] De maneira similar \[ \begin{array}{rcl} \mbox{Var}\big(\widehat{f}'(x)\big) & = & \mbox{Var}\left( \dfrac{\nu_0-\nu_{-1}}{nh^2}\right) = \dfrac{\mbox{Var}(\nu_0)-2\mbox{Cov}(\nu_0,\nu_{-1})+\mbox{Var}(\nu_{-1})}{n^2h^4}\\[0.8em] & = & \dfrac{p_0(1-p_0)-2(-p_0p_{-1}+p_{-1}(1-p_{-1}))}{nh^4} = \dfrac{2f(0)}{nh^3}+\cdots\, \cdot \end{array} \] Integrando sobre o intervalo, multiplica-se por \(h\) e somando sobre todos os intervalos, obtemos: \[ \mbox{AIV}=\dfrac{2}{nh^3}, \] desde que \(\displaystyle \int f(x)\mbox{d}x=1\).


Teorema 3.6: Suponha que \(f\) tenha uma segunda derivada absolutamente contínua que seja de quadrado integrável. Então, o \(\mbox{MISE}\) assintótico de \(\widehat{f}'(x)\) é \[ \mbox{AMISE}_{\widehat{f}'}(h) = \dfrac{2}{nh^3}+\dfrac{1}{12}h^2 R(f''); \] portanto, \[ \begin{array}{rcl} h^* & = & 6^{2/5} R(f'')^{-1/5}n^{-1/5}, \\[0.8em] \mbox{AMISE}^* & = & 5\times 6^{-6/5}R(f'')^{3/5}n^{-2/5}\cdot \end{array} \] Se \(f\sim N(\mu,\sigma^2)\), então \(h^*=2\times 3^{1/5}\pi^{1/10}\sigma\, n^{-1/5}\approx 2.8 \sigma\, n^{-1/5}\).


Demonstração Exercício.


Nas Figuras 3.25 e 3.26, são apresentadas porções de dez histogramas calculados a partir de amostras muito grandes de \(N(0,1)\), utilizando a largura de intervalo ideal para a densidade em comparação com a largura de intervalo ideal para a derivada da densidade.

Figura 3.25: Para 10 amostras \(N(0,1)\) de tamanho \(n = 107\), ampliações (deslocadas verticalmente) de porções dos 10 histogramas usando \(h^∗ = 0.01625\) na vizinhança de \(x = 0\), 1 e 3. Cada trecho do histograma mostra o intervalo de interesse e sete intervalos de cada lado.

Claramente, a largura de banda mais ampla fornece estimativas muito mais estáveis das inclinações na densidade verdadeira, embora ao custo de um viés maior. Onde a densidade é quase plana, muitos modos de amostra falsos são aparentes na Figura 3.26, enquanto todos os histogramas são monotonicamente decrescentes em uma vizinhança de \(x = 1\). Mesmo com tamanhos de amostra dessa magnitude, o problema da largura de banda deve sempre ser considerado potencialmente importante.

Figura 3.26: Para as mesmas 10 amostras \(N(0,1)\) de tamanho \(n = 10\), ampliações (deslocadas verticalmente) de porções dos 10 histogramas usando \(h^∗ = 0.1115\) na vizinhança de \(x = 0\), 1 e 3. Cada trecho do histograma mostra o intervalo de interesse e três intervalos de cada lado.

Na Figura 3.27, as estimativas da “derivada” do histograma são mostradas para toda a densidade \(N(0,1)\) para uma amostra de tamanho \(n = 10^5\), utilizando as larguras de banda ideais para a densidade e a derivada. Claramente, o ruído é bastante reduzido com o uso da largura de banda mais ampla. A largura de banda a ser escolhida depende da finalidade da análise.

Figura 3.27: Para uma amostra normal padrão com \(n = 10^5\) pontos, comparação das estimativas da “derivada” do histograma usando as larguras de banda ótimas de densidade e derivada.


3.5.4 Uma densidade de mistura bimodal útil


O estudo do desempenho de estimadores não paramétricos tem sido avaliado usando uma gama de densidades verdadeiras, por exemplo, as 15 misturas normais propostas por Marron and Wand (1992). No entanto, pode-se simplificar essa tarefa concentrando-se em apenas duas densidades: uma densidade normal padrão e uma mistura particular de duas normais dada por \[ \tag{3.78} f_M(x)=\dfrac{3}{4}\phi(x\, | \, 0,1)+\dfrac{1}{4}\phi(x\, | \, 0,1/3^2), \] que é uma versão simplificada de MW#14 (ver Figura 3.28). Essa densidade incorpora tanto a assimetria quanto diferentes graus de curvatura.

No entanto, o valor da densidade nas duas modas é praticamente idêntico. Mais adiante, essa densidade de teste se mostrará um desafio complexo em muitas situações.

Na Figura 3.28, são mostrados três histogramas para uma amostra de \(n = 1.000\) pontos da mistura. O primeiro utiliza a largura de banda ideal para o componente esquerdo e, portanto, o componente direito é suavizado em excesso.

Figura 3.28: Três histogramas de 1.000 pontos da mistura de dois componentes. As larguras de banda (da esquerda para a direita) são ótimas para o componente da esquerda, a mistura e o componente da direita, respectivamente.

O segundo utiliza a largura de banda ideal para a mistura (e, portanto, apresenta um desempenho ligeiramente pior para ambos). Finalmente, o terceiro histograma utiliza a largura de banda ideal para o componente direito (e, portanto, suaviza o componente esquerdo de forma insuficiente). Retornaremos a essa densidade de tempos em tempos.


3.6 Outros critérios de erro: \(L_1\), \(L_4\), \(L_6\), \(L_8\) e \(L_\infty\)


3.6.1 Histogramas \(L_1\) ótimos


Relembre os vários aspectos teoricamente atraentes do erro absoluto em vez do erro quadrático discutidos na Seção 2.3.2.2. A principal dificuldade com o erro absoluto é que ele não possui uma decomposição exata de variância-viés, embora Hall and Wand (1988a) tenham mostrado como estimar diretamente o erro total assintoticamente.

Hjort (1986) demonstrou que \[ \begin{array}{rcl} \mbox{E}\left( \int |\widehat{f}-f|\right) & \leq & \displaystyle \mbox{E}\left( \int |\widehat{f}-\mbox{E}\big(\widehat{f}\big)|\right) + \int |\mbox{E}\big(\widehat{f}\big)-f| \\[0.8em] & \approx & \displaystyle \sqrt{\dfrac{2}{\pi\, nh}}\int \sqrt{f} +\dfrac{1}{4}h\int |f '|\cdot \end{array} \] Portanto, minimizar esse limite superior assintótico para \(\mbox{MIAE}\) leva a \[ \tag{3.79} \begin{array}{rcl} h^* & = & \displaystyle 2\pi^{-1/3}\left(\int f^{1/2} \bigg/ \int |f'| \right)^{2/3}n^{-1/3} \\[0.8em] & = & \displaystyle 2.717 \sigma \, n^{-1/3} \quad \mbox{para dados } \; N(\mu,\sigma^2)\cdot \end{array} \]

No entanto, simulações numéricas com dados normais mostram que, na verdade, \(h^∗ ≈ 3.37\sigma \, n^{-1/3}\), que é 3.5% menor que a largura de banda \(L_2\) ideal. Minimizar os limites superiores do \(\mbox{MIAE}\) leva a larguras de banda que podem ser maiores ou menores que \(h^∗\).

Portanto, contrariamente à intuição, o erro absoluto nem sempre fornece bandas mais largas nas caudas. O que aconteceria com malhas \(L_1\) adaptativas ideais é desconhecido. Recentemente, Hall and Wand (1988a) examinaram as larguras de banda ideais para muitas densidades e descobriram que, com algumas densidades de cauda mais pesada, a largura de banda \(L_1\) pode ser maior que a largura de banda \(L_2\) ideal.


3.6.2 Outros critérios \(L_p\)


A análise de erros de histogramas para qualquer \(p\) par é direta usando momentos binomiais de ordem superior e procedendo como antes. Por exemplo, usando o erro \(L_4\), o erro médio integrado de quarta potência assintótico é (veja o Exercício 3.27) \[ \tag{3.80} \mbox{AMI4E}=\dfrac{3}{n^2h^2}\int f(x)^2\mbox{d}x+\dfrac{h}{2n}\int f'(x)f(x)\mbox{d}x+\dfrac{h^4}{80}\int f'(x)^4 \mbox{d}x\cdot \]

A Tabela 3.9 resume algumas larguras de banda ótimas conhecidas para dados \(N(0,1)\). A largura de banda ótima \(L_\infty\) é \(O(\log(n)n^{-1/3})\).

Tabela 3.9: Larguras de banda ideais para dados \(N(0,1)\) com diferentes critérios \(L_p\).

Novamente, qualquer critério de largura de banda fixa dará mais atenção às regiões onde a densidade é irregular; essa região não está necessariamente nas caudas da distribuição. Os coeficientes na tabela sugerem que \(L_4\) pode ter alguma vantagem especial, mas isso é um problema em aberto.


3.7 Exercícios


  • 1- Mostre como a regra da largura normal do intervalo pode ser modificada se \(f\) for assimétrica ou leptocurtótica, conforme discutido na introdução da Seção 3.3, usando outras densidades de referência. Examine o efeito da bimodalidade. Compare suas regras com as extensões da regra de Sturges propostas por Doane (1976).

  • 2- Realize uma análise dimensional para as quantidades \(f\),\(f'\),\(f''\), \(R(f)\), \(R(f')\) e \(R(f'')\). Verifique seus resultados calculando essas quantidades para o caso \(f = N(\mu,\sigma^2)\) rastreando o fator \(\sigma\).

  • 3- Uma aproximação para o erro de uma soma Riemanniana: \[ \left| \sum_{n=-\infty}^\infty g(nh)h -\int_{-\infty}^\infty g(x)\mbox{d}x \right| = \left|\sum_{n=-\infty}^{\infty} \int_{nh}^{nh+h} \Big(g(nh)-g(x) \Big)\mbox{d}x \right|\\ \leq \sum_{n=-\infty}^\infty \int_{nh}^{nh+h} \Big| g(nh)-g(x) \Big|\mbox{d}x \leq \sum_{n=-\infty}^\infty h V_g(nh,nh+h)\leq hV_g(\mathbb{R}^1), \] onde \(V_g(a,b)\) é a variação total de \(g\) em \([a,b]\) definida por \[ \sup \Big(\sum_{i=1}^n |g(x_i)-g(x_{i-1})| \Big), \] sobre todas as partições em \([a,b]\), incluindo \((a,b) = (-\infty,\infty)\). Conclua que se \(f'(\cdot)^2\) tem variação total finita, então o termo residual no viés (3.14) é \(O(h^3)\).

  • 4- Calcule a rugosidade de várias densidades paramétricas: Cauchy, \(t\)-Student, Beta e Lognormal. Para cada uma, calcule a largura de intervalo ótima. Expresse a largura de intervalo ótima em termos da variância, se existir. Compare os coeficientes nessas fórmulas com a regra normal na equação (3.16).

  • 5- Mostre que, quando \(h = h^*\) para o histograma, a contribuição dos termos \(\mbox{IV}\) e \(\mbox{ISB}\) para o \(\mbox{AMISE}\) é assintoticamente na proporção de 2:1.

  • 6- Verifique a equação (3.25).
    Dica: Substitua \(h = ch^∗\) na equação (3.23).

  • 7- Compare a sensibilidade do \(\mbox{AMISE}(ch^∗)\) na equação (3.25) para várias combinações de \(d\), \(p\) e \(r\).

  • 8- Demonstre as expressões exatas de \(\mbox{IV}\) e \(\mbox{ISB}\) na equação (3.26).

  • 9- Mostre que a \(\mbox{ISB}\) em um intervalo contendo a origem da densidade exponencial dupla, \[ f(x) = exp(-|x|)/2, \] é \(O(h^3)\); portanto, a descontinuidade na derivada de \(f\) não tem nenhum efeito assintótico na consistência. Compare a escolha \(t_0 = 0\) com a escolha \(t_0 = h/2\) como limite do intervalo. Formalmente, se \(f\) é a densidade exponencial negativa ordinária, \(R(f')\) é infinito devido ao salto em zero, integral do quadrado da função delta de Dirac em 0, mas \(R(f')\) está bem definido para a exponencial dupla.

  • 10- Considere o erro quadrático médio exato (\(\mbox{MISE}\)) de um histograma quando \(f = U(0,1)\).
      1. Se a malha \(t_k = kh\) for escolhida onde \(h = 1/m\), mostre que \[ \mbox{MISE}(m,n) = \dfrac{m-1}{n}\cdot \] Trivialmente, \(m^∗ = 1\) e \(h^∗ = 1\), e o erro é 0.

        1. Se \(t_k = (k+\dfrac{1}{2})h\) com \(h = 1/m\), mostre que \[ \mbox{MISE}(m,n)=\dfrac{1-2m+2m^2+n}{2mn}\cdot \] Conclua que \(h^*=\sqrt{2/(n+1)}\) e \[ \mbox{MISE}^*=\dfrac{\sqrt{2}\sqrt{n+1}-1}{n}\approx \sqrt{\dfrac{2}{n}}-\dfrac{1}{n}+\cdots, \] que é uma taxa mais lenta do que o histograma usual com densidade suave.

  • 11- Verifique as desigualdades para o \(\mbox{MISE}\) adaptativo na equação (3.32).

  • 12- Encontre as malhas adaptativas ótimas para uma densidade Beta assimétrica usando um programa de otimização numérica.

  • 13- Investigue o uso de malhas fixas e malhas percentis ao aplicar testes de hipótese de aderência qui-quadrado.

  • 14- Aplique o procedimento de sobresuavização ao conjunto de dados \(\mbox{LRL}\). Compare os resultados usando a amplitude e a variância como medidas de escala.

  • 15- Encontre uma regra de suavização excessiva baseada no intervalo interquartil como medida de escala.

  • 16- Use séries de Taylor e aproximações integrais de Riemann para verificar a equação (3.49).

  • 17- Avalie e verifique cuidadosamente as fórmulas na equação (3.53).
    Dica: \(\widehat{f}_{-i}(x_i)=(\nu_k-1)/\big((n-1)h\big)\).

  • 18- Quais são as larguras dos intervalos \(\mbox{UCV}\) e \(\mbox{BCV}\) para os conjuntos de dados \(\mbox{LRL}\) e de queda de neve (\(\mbox{snowfall}\))? Como os dados \(\mbox{LRL}\) já estão agrupados, tente várias larguras de intervalo que sejam múltiplos de 10 MeV e várias origens de intervalo.

  • 19- A regra \(\mbox{UCV}\) parece estimar algo entre \(\mbox{ISE}\) e \(\mbox{MISE}\), a quantidade intermediária \[ \displaystyle -2\int \widehat{f}(x)f(x)\mbox{d}x, \] sendo o foco. Usando a fórmula exata do \(\mbox{MISE}\) com dados normais e Lognormais, investigue o comportamento dos termos na aproximação \(\mbox{UCV}\).

  • 20- Considere uma amostra normal padrão de tamanho 1000. Aplique graus crescentes de arredondamento aos dados e compare as curvas \(\mbox{UCV}\) resultantes. Em que ponto o minimizador de \(\mbox{UCV}\) se torna \(h = 0\)?

  • 21- No \(\mbox{BCV}\), o viés na rugosidade foi reduzido por subtração. Alternativamente, a rugosidade poderia ter sido multiplicada por 3/4. Examine o efeito dessa ideia por meio de um exemplo. De onde se origina o fator 3/4?

  • 22- Desenvolva uma formulação \(\mbox{BCV}\) baseada em um estimador de diferença central de \(f'\) dado por \[ \dfrac{\dfrac{\nu_{k+1}}{nh}-\dfrac{\nu_k}{nh}}{2h}, \] ver Scott and Terrell (1987).

  • 23- Verifique as expressões nas equações (3.65) e (3.69) para as larguras dos intervalos e os erros quando a densidade for normal bivariada.

    • 24- Para dados normais bivariados, examine a ineficiência do uso de intervalos quadrados em relação aos intervalos retangulares. Examine combinações de situações em que a correlação é ou não é 0 e situações em que a razão das variâncias marginais é ou não é igual a 1.
  • 25- Verifique as expressões \(\mbox{AMISE}\) dadas nas equações (3.70) e (3.71) para as três malhas bivariadas regulares.
    Dica: Integre cuidadosamente o viés sobre a forma de cada intervalo e, em seguida, agregue usando a aproximação integral Riemanniana.

  • 26- Considere uma densidade de mistura normal de dois componentes, \[ f(x) = 0.5\times N(-\mu, 1) + 0.5\times N(\mu,1)\cdot \] Qual deve ser o valor de \(\mu\) para que a densidade apresente duas modas e picos?

  • 27- Calcule um dos critérios de erro com base em \(L_p\) para algum \(p > 2\) par.

3.8 Bibliografia


Bellman, R. E. 1961. Adaptive Control Processes. Princeton University Press, Princeton, NJ.
Box, G. E. P., and D. R. Cox. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society, B, no. 26: 211–43.
Carmichael, J.-P. 1976. “The Autoregressive Method: A Method for Approximating and Estimating Positive Functions.” PhD thesis, Ph.D. thesis, Statistical Science Division, SUNY, Buffalo, NY.
Carr, D. B., R. J. Littlefield, W. L. Nicholson, and J. S. Littlefield. 1987. “Scatterplot Matrix Techniques for Large n.” Journal of the American Statistical Association, no. 83: 596–610.
Day, N. E. 1969. “Estimating the Components of a Mixture of Normal Distributions.” Biometrika, no. 56: 463–74.
DIW. 1983. Das Sozio-Ökonomische Panel. Deutsches Institut für Wirtschafts-forschung, Berlin.
Doane, D. P. 1976. “Aesthetic Frequency Classifications.” American Statistician, no. 30: 181–83.
Donoho, D. L. 1988. “One-Sided Inference about Functionals of a Density.” Annals of Statistics, no. 16: 1390–420.
Duda, R. O., and P. E. Hart. 1973. Pattern Classification and Scene Analysis. John Wiley, New York.
Emerson, J. D., and D. C. Hoaglin. 1983. “Stem-and-Leaf Displays.” Edited by D. C. Hoaglin F. Mosteller and J. W. Tukey. John Wiley, New York.
Epanechnikov, V. K. 1969. “Non-Parametric Estimation of a Multivariate Probability Density.” Theory of Probability and Its Applications, no. 14: 153–58.
Everitt, B. S., and D. J. Hand. 1981. Finite Mixture Distributions. Chapman; Hall, London.
Freedman, D., and P. Diaconis. 1981. “On the Histogram as a Density Estimator: L2 Theory.” Zeitschrift Für Wahrscheinlichkeitstheorie Und Verwandte Gebiete, no. 57: 453–76.
Good, I. J., and R. A. Gaskins. 1980. “Density Estimation and Bump-Hunting by the Penalized Likelihood Method Exemplified by the Scattering and Meteorite Data (with Discussion).” Journal of the American Statistical Association, no. 75: 42–73.
Hall, P., and J. S. Marron. 1987a. “Estimation of Integrated Squared Density Derivatives.” Statistics and Probability Letters, no. 6: 109–15.
———. 1987b. “Extent to Which Least-Squares Cross-Validation Minimises Integrated Square Error in Nonparametric Density Estimation.” Probability Theory and Related Fields, no. 74: 567–81.
———. 1987c. “On the Amount of Noise Inherent in Bandwidth Selection for a Kernel Density Estimator.” Annals of Statistics, no. 15: 163–81.
Hall, P., and M. P. Wand. 1988a. “Minimizing L1 Distance in Nonparametric Density Estimation.” Journal of the Multivariate Analysis, no. 26 (1988a): 59–88.
Hathaway, R. J. 1985. “A Constrained Formulation of Maximum-Likelihood Estimation for Normal Mixture Distributions.” Annals of the Statistics, no. 13: 795–800.
Hjort, N. L. 1986. “On Frequency Polygons and Averaged Shifted Histograms in Higher Dimensions.” Stanford University.
Hodges, J. L., and E. L. Lehmann. 1956. “The Efficiency of Some Nonparametric Competitors of the t-Test.” Annals of the Mathematical Statistics, no. 27: 324–35.
Hüsemann, J. A., and G. R. Terrell. 1991. “Optimal Parameter Choice for Error Minimization in Bivariate Histograms.” Journal of Multivariate Analysis, no. 37: 85–103.
Izenman, A. J., and C. J. Sommer. 1988. “Philatelic Mixtures and Multimodal Densities.” Journal of the American Statistical Association, no. 83: 941–53.
Kogure, A. 1987. “Asymptotically Optimal Cells for a Histogram.” Annals of Statistics, no. 15: 1023–30.
Marron, J. S., and M. P. Wand. 1992. “Exact Mean Integrated Squared Error.” Annals of Statistics, no. 20: 712–36.
Parzen, E. 1979. “Nonparametric Statistical Data Modeling.” Journal of the American Statistical Association, no. 74: 105–31.
Redner, R. A., and H. F. Walker. 1984. “Mixture Densities, Maximum Likelihood and the EM Algorithm.” SIAM Review, no. 26: 195–202.
Rudemo, M. 1982. “Empirical Choice of Histograms and Kernel Density Estimators.” Scandinavian Journal of Statistics, no. 9: 65–78.
Scott, D. W. 1979. “On Optimal and Data-Based Histograms.” Biometrika, no. 66: 605–10.
———. 1988. “A Note on Choice of Bivariate Histogram Bin Shape.” Journal of Official Statistics, 47–51.
Scott, D. W., and H.-P. Schmidt. 1988. “Calibrating Histograms with Applications to Economic Data.” Empirical Economics, no. 13: 155–68.
———. 1989. “Calibrating Histograms with Applications to Economic Data.” Edited by A. Ullah. Physica-Verlag, Heidelberg.
Scott, D. W., and G. R. Terrell. 1987. “Biased and Unbiased Cross-Validation in Density Estimation.” Journal of the American Statistical Association, no. 82: 1131–46.
Scott, D. W., and J. R. Thompson. 1983. “Probability Density Estimation in Higher Dimensions.” In Proceedings of the Fifteenth Interface of Computer Science and Statistics, edited by J. E. Gentle, 173–79. North-Holland, Amsterdam.
Terrell, G. R. 1990. “The Maximal Smoothing Principle in Density Estimation.” Journal of the Amererican Statistical Association, no. 85: 470–77.
Terrell, G. R., and D. W. Scott. 1983. “Variable Window Density Estimates.” Technical report presented at ASA meetings in Toronto.
———. 1985. “Oversmoothed Nonparametric Density Estimates.” Journal of the American Statistical Association, no. 80.
———. 1992. “Variable Kernel Density Estimation.” Annals of Statistics, no. 20: 1236–65.
Tukey, J. W. 1977. Exploratory Data Analysis. Addison-Wesley, Reading, MA.
Van Ryzin, J. 1973. “A Histogram Method of Density Estimation.” Communications in Statistics, no. 2: 493–506.
Wand, M. P., and M. C. Jones. 1991. “Comparison of Smoothing Parameterizations in Bivariate Density Estimation.” Journal of the American Statistical Association, no. 88: 520–28.
Wegman, E. J. 1970. “Maximum Likelihood Estimation of a Unimodal Density Function.” Annals of Statistics, no. 41: 457–71.
Wharton, S. W. 1983. “A Generalized Histogram Clustering Scheme for Multidimensional Image Data.” Pattern Recognition, no. 16: 193–99.