Capítulo 6

Estimadores kernel de densidade


É notável que o histograma tenha permanecido como o único estimador de densidade não paramétrico até a década de 1950, quando progressos substanciais e simultâneos foram feitos na estimação de densidade e na estimação de densidade espectral. Em um relatório técnico pouco conhecido, Fix and Hodges (1951) introduziram o algoritmo básico da estimação de densidade não paramétrica.

Eles abordaram o problema da discriminação estatística quando a forma paramétrica da densidade de amostragem não era conhecida. Felizmente, este artigo foi reimpresso com comentários por Silverman and Jones (1989). Durante a década seguinte, diversos algoritmos gerais e modos alternativos de análise teórica foram introduzidos por Rosenblatt (1956), Parzen (1962) e Cencov (1962).

Seguiu-se uma segunda onda de artigos importantes e principalmente teóricos de Watson and Leadbetter (1963), Loftsgaarden and Quesenberry (1965), Schwartz (1967), Epanechnikov (1969), Tarter and Kronmal (1970) e Wahba (1971). A generalização multivariada natural foi introduzida por Cacoullos (1966). Finalmente, na década de 1970, surgiram os primeiros artigos focados na aplicação prática desses métodos: Scott et al. (1978) e Silverman (1978b). Essas e outras aplicações multivariadas posteriores aguardavam a revolução da computação.

O estimador de kernel básico pode ser escrito de forma compacta como \[ \tag{6.1} \widehat{f}(x)=\dfrac{1}{nh}\sum_{i=1}^n K\left( \dfrac{x-x_i}{h}\right)=\dfrac{1}{n}\sum_{i=1}^n K_h (x-x_i), \] onde \(K_h(t)=K(t/h)/h\), que é uma notação introduzida por Rosenblatt (1956).

O estimador kernel pode ser motivado não apenas como o caso limite do histograma deslocado médio como em (5.14), mas também por outras técnicas demonstradas na Seção 6.1.

De fato, praticamente todos os algoritmos não paramétricos são assintoticamente métodos kernel, um fato demonstrado empiricamente por Walter and Blum (1979a) e provado rigorosamente por G. R. Terrell and Scott (1992a). Woodroofe (1970) chamou a classe geral de “sequências delta”.


  • 6.1 Motivação para estimadores kernel
  • 6.2 Propriedades teóricas: caso univariado
  • 6.3 Propriedades teóricas: caso multivariado
  • 6.4 Generalidade do método kernel
  • 6.5 Validação cruzada
  • 6.6 Suavização adaptativa
  • 6.7 Aspectos computacionais

    • 6.1 Motivação para estimadores kernel


    Do ponto de vista de um estatístico ou instrutor, a média de histogramas deslocados parece a motivação mais natural para estimadores de kernel. No entanto, seguir outros pontos de partida em análise numérica, séries temporais e processamento de sinais proporciona uma compreensão mais profunda dos métodos de kernel.

    Ao tentar compreender um ponto teórico ou prático específico relativo a um estimador não paramétrico, nem todas as abordagens são igualmente poderosas. Por exemplo, a análise de Fourier fornece ferramentas sofisticadas para fins teóricos. A relação entre viés e variância pode ser reformulada em termos de filtros passa-baixo e passa-alto no processamento de sinais. Cada um descreve a mesma entidade, mas com matemática diferente.


    6.1.1 Análise numérica e diferenças finitas


    O estimador kernel surgiu como uma aproximação numérica da derivada da função de distribuição acumulada (Rosenblatt 1956). A função de densidade de probabilidade empírica, definida na equação (2.2) como a derivada formal da função de distribuição acumulada empírica \(F_n(x)\), é uma soma de funções delta de Dirac, o que a torna inútil como estimador de uma função de densidade suave.

    Considere, no entanto, uma aproximação por diferenças finitas unilateral para a derivada de \(F_n(\cdot)\): \[ \tag{6.2} \begin{array}{rcl} \widehat{f}(x) & = & \displaystyle \dfrac{F_n(x)-F_n(x-h)}{h}\\[0.8em] & = & \displaystyle \dfrac{1}{nh}\sum_{i=1}^n \pmb{I}_{[x-h,x)}(x_i) =\dfrac{1}{nh} \sum_{i=1}^n \pmb{I}_{(0,1]}\left( \dfrac{x-x_i}{h}\right), \end{array} \] que, a partir da equação (6.1), é claramente um estimador de kernel com \(K = U(0,1]\).

    Como \(\mbox{E}\big(F_n(x)\big) = F(x)\) para todo \(x\), então, com a série de Taylor \[ F(x-h)=F(x)-h f(x)+\dfrac{1}{2}h^2 f'(x)-\dfrac{1}{6}h^3 f''(x)+ \cdots, \] o viés é calculado como \[ \mbox{Viés}\big( \widehat{f}(x)\big)=\mbox{E}\big(\widehat{f}(x)-f(x) \big)=-\dfrac{1}{2}hf'(x)+O(h^2)\cdot \] Assim, o viés quadrático integrado é \(h^2 R(f')/4\), que é comparável à ordem do viés quadrático integrado (\(\mbox{ISB}\)) do histograma no Teorema 3.1, em vez do \(\mbox{ISB}\) \(O(h^4)\) do polígono de frequência (FP).

    Além disso, o \(\mbox{ISB}\) de (6.2) é três vezes maior que o \(\mbox{ISB}\) do histograma. As integrais das variâncias são idênticas; veja o Exercício 1. Portanto, o estimador de kernel unilateral (6.2) é inferior a um histograma.

    Sem comentários, Rosenblatt propôs um estimador bilateral ou de diferença central de \(f\): \[ \tag{6.3} \widehat{f}(x)\dfrac{F_n(x+h/2)-F_n(x-h/2)}{h}\cdot \] O viés de (6.3) resulta em \(h^2 f''(x)/24\), ver Exercício 2. Assim, o viés ao quadrado é \(O(h^4)\), o mesmo que o do FP. O kernel correspondente é \(K = U(-0.5; 0.5)\).

    Lembre-se de que o histograma colocou um bloco retangular no compartimento onde cada ponto de dados se encontrava. O estimador unilateral (6.2) coloca a borda esquerda de um bloco retangular em cada ponto de dados, enquanto o estimador bilateral (6.3) coloca o centro de um bloco retangular em cada ponto de dados (ver Tarter and Kronmal 1976).

    A Figura 6.1 exibe duas estimativas de diferença central dos dados de queda de neve de Buffalo. Comparadas com o FP, essas estimativas são graficamente inferiores, pois a estimativa contém \(2n\) saltos que nem sequer são igualmente espaçados. No entanto, a maioria dos critérios, como o \(\mbox{MISE}\), não é particularmente sensível a esse comportamento local ruidoso.

    Figura 6.1: Estimativas de diferença central dos dados de queda de neve em Buffalo.

    Os códigos de otimização quase-Newton rotineiramente utilizam estimativas numéricas de diferença central de derivadas. Alguns códigos usam aproximações de “ordem ainda maior” para a primeira derivada (ver Seção 6.2.3.1).


    6.1.2 Suavização por convolução


    Um engenheiro eletricista que se depara com uma função ruidosa recorrerá a uma variedade de filtros de convolução para encontrar um que suavize os componentes indesejados de alta frequência. A operação de convolução “∗”, que é definida como \[ (f * \omega)(x)=\int_{-\infty}^\infty f(t)\omega(x-t)\mbox{d}t, \] substitui o valor de uma função \(f(x)\), por uma média ponderada local dos valores da função, de acordo com uma função de ponderação \(\omega(\cdot)\) que geralmente é simétrica e concentrada em torno de 0. Os estatísticos também utilizam a operação de média para reduzir a variância.

    Portanto, a função de densidade de probabilidade empírica (2.2), que é muito ruidosa, pode ser filtrada por convolução, com o resultado de que \[ \tag{6.4} \begin{array}{rcl} \left(\dfrac{\mbox{d}F_n}{\mbox{d}x}\right) * \omega & = & \displaystyle \int_{-\infty}^\infty \left(\dfrac{1}{n}\sum_{i=1}^n \delta(t-x_i) \right)\omega(x-t)\mbox{d}t \\[0.8em] & = & \displaystyle \dfrac{1}{n}\sum_{i=1}^n \left(\int_{-\infty}^\infty \delta(t-x_i)\, \omega(x-t)\mbox{d}t \right) = \dfrac{1}{n}\sum_{i=1}^n \omega(x-x_ì), \end{array} \] que é a segunda forma de kernel dada em (6.1), mas sem o parâmetro de suavização \(h\) aparecendo explicitamente como, por exemplo, \(\omega(t) = K_h(t)\).

    Em geral, a forma e a extensão da função de ponderação do filtro de convolução \(\omega\) dependerão do tamanho da amostra. O estimador de kernel (6.1) usa uma única “forma” para todos os tamanhos de amostra, e a largura do kernel é controlada explicitamente pelo parâmetro de suavização \(h\). A literatura sobre projeto de filtros frequentemente usa terminologia diferente. Por exemplo, a largura do filtro \(\omega\) pode ser controlada pelo ponto de meia potência, onde o filtro atinge metade do seu valor na origem.


    6.1.3 Aproximações por séries ortogonais


    A introdução heurística ao alisamento por convolução pode ser formalizada por um argumento de aproximação por séries ortogonais. Por simplicidade, suponha que a função densidade \(f\) seja periódica no intervalo [0,1], de modo que a base de séries de Fourier ordinárias \(\phi_\nu (t) = \exp\big(2\pi i\nu t\big)\), seja apropriada.

    Toda função, mesmo funções ruidosas, pode ser expressa em termos dessas funções de base como \[ \tag{6.5} f(x)=\sum_{\nu=-\infty}^\infty f_\nu \phi_\nu (x), \qquad \mbox{onde} \qquad f_\nu=<f,\phi_\nu>=\int_0^1 f(x)\phi_\nu^*(x)\mbox{d}x\cdot \] As funções de base são ortonormais, isto é, \[ \int \phi_\nu^*(x)\phi_\mu(x)\mbox{d}x = \delta_{\nu\mu}, \] onde \(\delta_{\nu\mu}\) é a função delta de Kronecker e \(\phi^∗\) denota o complexo conjugado.

    Como \(f\) é uma função de densidade, o coeficiente \(f_\nu\) na equação (6.5) pode ser expresso em termos estatísticos como \[ \tag{6.6} f_\nu=\mbox{E}\big(\phi_\nu^*(X) \big); \qquad \mbox{então} \qquad \widehat{f}_\nu =\dfrac{1}{n}\sum_{\ell=1}^n \phi_\nu^*(x_\ell), \] é um estimador não viesado e consistente do coeficiente de Fourier \(f_\nu\). Observe que a soma é sobre \(\ell\) para evitar confusão com \(i = \sqrt{-1}\).

    Como um exemplo extremo, considere os coeficientes de Fourier da função de densidade de probabilidade empírica (2.2): \[ f_\nu =\int_0^1 \left( \dfrac{1}{n}\sum_{\ell=1}^n \delta(x-x_\ell)\right)\phi_\nu^*(x)\mbox{d}x=\dfrac{1}{n}\sum_{\ell=1}^n \phi_\nu^*(x_\ell)=\widehat{f}_\nu, \] onde é \(\widehat{f}_\nu\) definido em (6.6).

    Como \(\widehat{f}_\nu\) e \(f_\nu\) para a função densidade de probabilidade empírica são idênticos, o seguinte é formalmente verdadeiro para qualquer amostra \(\{x_\ell\}\): \[ \tag{6.7} \sum_{\nu=-\infty}^\infty \widehat{f}_\nu \phi_\nu(x)=\dfrac{1}{n}\sum_{\ell=1}^n \delta(x-x_\ell), \] que é a função de densidade de probabilidade empírica.

    Cencov (1962), Kronmal and Tarter (1968) e Watson (1969) sugeriram suavizar a função de densidade empírica incluindo apenas alguns termos selecionados da equação (6.7).

    Excluir termos da forma \(|\nu|> k\) corresponde ao que os engenheiros chamam de pesos de filtro “boxcar” \[ \tag{6.8} \omega_\nu(k)=\left\{ \begin{array}{rcl} 1, & \mbox{caso}, & |\nu|> k \\[0.8em] 0, & \mbox{caso} & |\nu|\leq k \end{array}\right.\cdot \]

    Como a transformada de Fourier da função boxcar é a função \(\mbox{sinc}\), \(\sin(x)/(\pi x)\), a estimativa será grosseira e apresentará “vazamento”; ou seja, pontos amostrais relativamente distantes de um ponto \(x\) influenciarão \(\widehat{f}(x)\).

    Wahba (1977) sugere a aplicação de uma janela de suavização a esta série, o que proporciona um ajuste mais preciso da estimativa resultante. Ela introduz dois parâmetros, \(\lambda\) e \(\rho\), que controlam a forma e a extensão da janela de suavização: \[ \tag{6.9} \omega_\nu(\lambda,\rho)=\dfrac{1}{1+\lambda(2\pi \nu)^{2\rho}}, \qquad \mbox{para} \qquad |\nu|\leq n/2\cdot \] Ambas as formas da estimativa de Fourier ponderada podem ser escritas explicitamente como \[ \tag{6.10} \widehat{f}_\nu(x)=\sum_\nu \omega_\nu \left(\dfrac{1}{n}\sum_{\ell=1}^n \phi_\nu^* (x_\ell) \right)\phi_\nu(x)=\dfrac{1}{n}\sum_{\ell=1}^n \left(\sum_\nu \omega_\nu \phi_\nu^*(x_\ell)\phi_\nu(x) \right), \] onde a ordem das somas foi trocada.

    Com a base de Fourier, o estimador de série ortogonal (6.10) é igual a \[ \widehat{f}(x)=\dfrac{1}{n}\sum_{\ell=1}^n \left(\sum_\nu \omega_\nu e^{2\pi \ i\nu (x-x_\ell)} \right)\cdot \] Este estimador encontra-se agora na forma de convolução (6.4) de um estimador de kernel fixo, com o filtro (ou kernel) definido pela quantidade entre parênteses.

    Alguns exemplos dessas funções kernel são mostrados na Figura 6.2. Os kernels equivalentes de Wahba são mais suaves e apresentam menor vazamento. Os parâmetros \(\lambda\) e \(\rho\) podem ser interpretados como correspondentes ao parâmetro de suavização e à ordem da aproximação por diferenças finitas, respectivamente.

    Figura 6.2: Exemplos de kernels equivalentes para estimadores de séries ortogonais. Os quatro kernels de Wahba (linha inferior) foram selecionados para corresponder à altura do pico dos kernels de Kronmal-Tarter-Watson correspondentes (linha superior). Os kernels de Kronmal-Tarter-Watson são independentes do tamanho da amostra; os exemplos de Wahba são para \(n = 16\).


    6.2 Propriedades teóricas: caso univariado


    6.2.1 Análise \(\mbox{MISE}\)


    A análise estatística de estimadores kernel é muito mais simples do que a de histogramas, visto que o estimador de kernel (6.1) é a média aritmética de \(n\) variáveis aleatórias independentes e identicamente distribuídas \[ K_h(x,X_i)=\dfrac{1}{h}K\left( \dfrac{x-X_i}{h}\right)\cdot \] Portanto, \[ \tag{6.11} \mbox{E}\big(\widehat{f}(x) \big)=\mbox{E}\big( K_h(x,X)\big) \qquad \mbox{e} \qquad \mbox{Var}\big(\widehat{f}(x) \big)=\dfrac{1}{n}\mbox{Var}\big( K_h(x,X)\big)\cdot \] A esperança é igual a \[ \tag{6.12} \mbox{E}\big( K_h(x,X)\big) = \displaystyle \int \dfrac{1}{h}K\left(\dfrac{x-t}{h} \right)f(t)\mbox{d}t = \int K(\omega) f(x-h\omega)\mbox{d}\omega \] que pode ser escrita como \[ \tag{6.13} \mbox{E}\big( K_h(x,X)\big) = \displaystyle f(x)\int K(\omega)\mbox{d}\omega - hf'(x)\int \omega K(\omega)\mbox{d}\omega +\dfrac{1}{2}h^2 f''(x)\int \omega^2 K(\omega)\mbox{d}\omega+\cdots , \] e a variância é dada por \[ \tag{6.14} \mbox{Var}\big( K_h(x,X)\big) = \displaystyle \mbox{E}\left(\dfrac{1}{h}K\left(\dfrac{x-X}{h} \right) \right)^2-\left(\mbox{E}\left(\dfrac{1}{h}K\left(\dfrac{x-X}{h} \right) \right)\right)^2\cdot \] O segundo termo em (6.14) foi calculado em (6.13) e é aproximadamente igual a \[ \left(f(x)\int K(\omega)\mbox{d}\omega +\cdots \right)^2, \] enquanto o primeiro termo pode ser aproximado por \[ \tag{6.15} \int \dfrac{1}{h^2}K\left(\dfrac{x-X}{h}\right)f(t)\mbox{d}t=\int \dfrac{1}{h}K(\omega)^2 f(x-h\omega)\mbox{d}\omega\approx \dfrac{f(x)R(K)}{h}\cdot \]

    Da equação (6.13), se o núcleo \(K\) satisfizer \[ \int K(\omega)\mbox{d}\omega=1, \qquad \int \omega K(\omega)\mbox{d}\omega=0 \qquad \mbox{e} \qquad \int \omega^2 K(\omega)\mbox{d}\omega=\sigma_K^2>0, \] então a esperança de \(\widehat{f}(x)\) será igual a \(f(x)\) até a ordem \(O(h^2)\). De fato, \[ \tag{6.16} \mbox{Viés}(x)=\dfrac{1}{2}\sigma_K^2 h^2 f''(x)+O(h^4) \qquad \mbox{e} \qquad \mbox{ISB}=\dfrac{1}{4}\sigma_K^4 h^4 R(f'')+O(h^6)\cdot \] Da mesma forma, a partir de (6.14), (6.15) e (6.13), \[ \tag{6.17} \begin{array}{rcl} \mbox{Var}(x) & = & \displaystyle \dfrac{f(x)R(K)}{nh}-\dfrac{f(x)^2}{n}+O(h/n),\\[0.8em] \mbox{IV} & = & \displaystyle \dfrac{R(K)}{nh}-\dfrac{R(f)}{n}+\cdots\,\cdot \end{array} \] Esses resultados estão resumidos no seguinte teorema.


    Teorema 6.1: Para um estimador de densidade de kernel univariado não negativo, \[ \tag{6.18} \begin{array}{rcl} \mbox{AMISE} & = & \displaystyle \dfrac{R(K)}{nh}+\dfrac{1}{4}\sigma_K^4 h^4 R(f''), \\[0.8em] h^* & = & \displaystyle \left(\dfrac{R(K)}{\sigma_K^4 R(f'')} \right)^{1/5} n^{-1/5}\\[0.8em] \mbox{AMISE}^* & = & \displaystyle \dfrac{5}{4}\Big( \sigma_K R(K)\Big)^{4/5} R(f'')^{1/5} n^{-4/5}\cdot \end{array} \]


    Demosntração. As condições sob as quais o teorema é válido foram exploradas por muitos autores, incluindo Parzen (1962).


    Um conjunto simples de condições é que o núcleo \(K\) seja uma função de densidade de probabilidade contínua com suporte finito, \(K\in L_2\), \(\mu_K = 0\), \(0 <\sigma_K^2 <\infty\), que \(f\) seja absolutamente contínua e \(f'''\in L_2\) (Scott 1985).

    É possível verificar que a razão entre a integral da variância assintótica (\(\mbox{AIV}\)) e o viés quadrático integrado assintótico (\(\mbox{AISB}\)) no erro quadrático integrado médio assintótico (\(\mbox{AMISE}^*\)) é de 4:1. Ou seja, o \(\mbox{ISB}\) compreende apenas 20% do \(\mbox{AMISE}\).

    A semelhança com os resultados \(\mbox{FP}\) no Teorema 4.1 é clara. Se \(K\) for um triângulo isósceles, então os resultados no Teorema 6.1 correspondem aos do histograma médio deslocado (\(\mbox{ASH}\)) ingênuo com \(m =\infty\) na equação (5.8) (ver Exercício 3).

    Como \[ R\big(\phi'' (x\, | \, 0.\sigma^2)\big) = \dfrac{3}{8\sqrt{\pi}\sigma^5}, \] a largura de banda da regra de referência normal com um kernel normal é \[ \tag{6.19} \mbox{largura de banda da regra de referência normal} = (4/3)^{1/5}\sigma \, n^{-1/5}\approx 1.06 \, \widehat{\sigma}\, n^{-1/5}\cdot \]


    6.2.2 Estimação de derivadas


    Ocasionalmente, surge a necessidade de estimar as derivadas da função de densidade; por exemplo, ao procurar modas e picos. As derivadas de uma estimativa de kernel comum comportam-se de forma consistente se o kernel for suficientemente diferenciável e se forem selecionadas larguras de banda maiores. Parâmetros de suavização maiores são necessários, pois a derivada da função estimada é mais ruidosa do que a própria função estimada.

    Considere como estimador da \(r\)-ésima derivada de \(f\) a \(r\)-ésima derivada da estimativa de kernel: \[ \tag{6.20} \widehat{f}^{(r)}(x)= \dfrac{\mbox{d}^r}{\mbox{d}x^r} \dfrac{1}{nh}\sum_{i=1}^n K\left( \dfrac{x-x_i}{h}\right)=\dfrac{1}{nh^{r+1}} \sum_{i=1}^n K^{(r)}\left( \dfrac{x-x_i}{h}\right)\cdot \] Um cálculo semelhante ao que levou à equação (6.15) mostra que \[ \mbox{Var}\big(\widehat{f}^{(r)}(x)\big)\approx \dfrac{n}{\big(nh^{r+1} \big)^2} \mbox{E}\left( K^{(r)}\left(\dfrac{x-X}{h} \right)^2\right)\approx \dfrac{f(x)R(K^{(r)})}{nh^{2r+1}}, \] portanto, a integral da variância assintótica de \(\widehat{f}^{(r)}\) é \(R(K^{(r)})/(nh^{2r+1})\).

    Após uma expansão semelhante à de (6.13) para encontrar o viés, a esperança do estimador da primeira derivada é \[ \mbox{E}\big(\widehat{f}^{(1)}(x)\big)=\dfrac{1}{h}\left( f_x \int K' -hf_x' \int\omega K'+\dfrac{h^2}{2}f_x'' \int \omega^2 K' -\dfrac{h^3}{6} f_x'''\int \omega^3 K'+\cdots\right), \] onde \(f_x^{(r)}=f^{(r)}(x)\).

    Supondo que \(K\) seja simétrico, \(\displaystyle \int \omega^r K' = 0\) para \(r\) par. Integrando por partes, \(\displaystyle \int \omega K'= −1\) e \(\displaystyle \int \omega^3 K' = -3\sigma_K^2\). Portanto, o viés pontual é de ordem \(h^2\) e envolve a terceira derivada de \(f\). Um teorema geral é obtido (veja o Exercício 6).


    Teorema 6.2: Com base em um estimador de densidade kernel univariado não negativo \(\widehat{f}\), \[ \tag{6.21} \mbox{AMISE}\big(\widehat{f}^{(r)}(x)\big)=\dfrac{R(K^{(r)})}{nh^{2r+1}}+\dfrac{1}{4} h^4 \sigma_K^4 R(f^{(r+2)}), \] \[ h_r^* = \left(\dfrac{(2r+1)R(K^{(r)})}{\sigma_K^4 R(f^{(r+2)})} \right)^{1/(2r+5)} n^{-1/(2r+5)} \] e \[ \mbox{AMISE}^* \big(\widehat{f}^{(r)}(x)\big)=\dfrac{2r+5}{4}R(K^{(r)})^{4/(2r+5)}+\left(\sigma_K^4 R(f^{r+2}/(2r+1)) \right)^{(2r+1)/(2r+5)} n^{-4/(2r+5)}\cdot \]


    Demonstração. Scott (2015).


    Embora a ordem do termo de viés permaneça \(O(h^4)\), cada ordem de derivada adicional introduz duas potências extras de \(h\) na variância. Os parâmetros de suavização ótimos \(h^*\) para a primeira e a segunda derivadas são de ordem \(O(n^{-1/7})\) e \(O(n^{-1/9})\), respectivamente, enquanto o \(\mbox{AMISE}^*\) é de ordem \(O(n^{-4/7})\) e \(O(n^{-4/9})\).

    Se a taxa de densidade ótima \(h^* = O(n^{-1/5})\) for usada na estimativa da segunda derivada, a \(\mbox{IV}\) assintótica no \(\mbox{AMISE}\) não desaparece, uma vez que \(nh^5 = O(1)\). A estimativa de uma derivada adicional é mais difícil do que estimar uma dimensão adicional. Por exemplo, a taxa \(\mbox{AMISE}\) ótima para a segunda derivada é \(O(n^{-4/9})\), que é a mesma taxa (mais lenta) que a do \(\mbox{AMISE}\) ótimo de um estimador de densidade de polígono de frequência multivariado de 5 dimensões.


    6.2.3 Escolha do kernel


    Grande parte da primeira década de trabalho teórico concentrou-se em vários aspectos das propriedades de estimação relacionadas às características de um kernel ou núcleo. Dentro de uma classe particular de kernel, por exemplo, a ordem do seu primeiro momento não nulo, a qualidade de uma estimativa de densidade é agora amplamente reconhecida como sendo determinada principalmente pela escolha do parâmetro de suavização e apenas de forma secundária pela escolha do kernel, como ficará evidente na Tabela 6.2.

    Assim, o tópico poderia ser menos enfatizado. No entanto, houve um recente aumento de pesquisas úteis sobre o projeto de kernels em situações especiais. Embora muitos riscos potenciais afetem o usuário da estimação de densidade, por exemplo, subestimar a suavidade da densidade desconhecida, a especificação das propriedades desejadas para o kernel está inteiramente à disposição do pesquisador, que deve ter um bom entendimento dos resultados a seguir.


    6.2.3.1 Kernels de ordem superior


    Bartlett (1963) considerou a possibilidade de escolher cuidadosamente o núcleo para reduzir ainda mais a contribuição do viés ao \(\mbox{MISE}\). Se a exigência de que a estimativa do núcleo seja uma densidade verdadeira for relaxada, então é possível obter uma melhoria significativa no \(\mbox{MISE}\).

    Suponha que um núcleo de ordem \(p\) seja escolhido de forma que \[ \tag{6.22} \int K =1, \qquad \int \omega^i K=0, \quad i=1,\cdots,p-1; \qquad \mbox{e} \qquad \int \omega^p K\neq 0, \] então, continuando a expansão na equação (6.13), o viés pontual do kernel torna-se \[ \mbox{Viés}\big(\widehat{f}(x)\big)=\dfrac{1}{p!}h^p \mu_p \, f^{(p)}(x)+\cdots, \] considerando \(\mu_i = \displaystyle \int \omega^i K(\omega)\mbox{d}\omega\).

    Como as fórmulas para as variâncias pontuais e integradas permanecem inalteradas, o seguinte teorema pode ser provado.


    Teorema 6.3: Supondo que \(f\) seja suficientemente diferenciável e que o núcleo \(K\) seja de ordem \(p\), \[ \mbox{AMISE}(h)=\dfrac{R(K)}{nh}+\dfrac{1}{(p!)^2} \mu_p^2 \, R(f^{(r)}) h^{2p}, \] \[ \tag{6.23} h^* = \left( \dfrac{(p!)^2 R(K)}{2p\mu_p^2 \, R(f^{(p)})}\right)^{1/(2p+1)} n^{-1/(2p+1)} \] e \[ \mbox{AMISE}^*=\dfrac{2p+1}{2p}\left(2p \, \mu_p^2 \, R(K)^{2p} R(f^{(p)})/(p!)^2 \right)^{2/(2p+1)} n^{-2p/(2p+1)}\cdot \]


    Demonstração. Scott (2015).


    Assintoticamente, este resultado indica que é possível aproximar-se da taxa paramétrica usual de ordem \(O(n^{-1})\) para o \(\mbox{AMISE}\). No entanto, a largura das bandas ótimas aumenta à medida que a ordem do kernel \(p\) aumenta, sugerindo que grande parte do benefício pode ser bastante assintótica para valores grandes de \(p\).

    A Tabela 6.1 mostra alguns kernels de ordem superior, juntamente com o \(\mbox{AMISE}\) ótimo para o caso de dados normais padrão. Esses kernels são mostrados na Figura 6.3. Selecionar kernels representativos de ordem superior é difícil, pelas razões apresentadas na Seção 6.2.3.2. Apenas valores pares de \(p\) são considerados, porque todos os “momentos” ímpares de kernels simétricos se anulam.

    Tabela 6.1: Alguns núcleos polinomiais simples de ordem superior.

    Cada aumento de 2 em \(p\) adiciona outra restrição de momento (par) na equação (6.22). Se os kernels forem polinômios, o grau do polinômio também deve aumentar para que haja graus de liberdade suficientes para satisfazer as restrições. Os kernels na Tabela 6.1 começam com o chamado kernel de Epanechnikov e são os únicos kernels polinomiais contínuos de grau \(p\) que satisfazem as restrições e têm seu suporte no intervalo [-1,1].

    Os gráficos do \(\mbox{AMISE}\) para dados normais na Figura 6.3 sugerem que os kernels de ordem superior requerem vários milhares de pontos de dados antes que um ganho substancial possa ser obtido. Para dados mais irregulares, os ganhos são ainda mais assintóticos. A melhoria possibilitada pela utilização de um kernel de ordem superior não é simplesmente um fator multiplicativo constante, mas sim uma mudança exponencial na ordem de convergência de \(n\).

    Figura 6.3: Exemplos de kernels de ordem superior que são polinômios de ordem inferior. O painel da direita mostra as curvas \(\mbox{AMISE}^*\) da distribuição \(N(0,1)\) correspondentes, em escala log-log.

    Obviamente, para amostras pequenas, a diferença entre \(\mbox{MISE}\) e \(\mbox{AMISE}\) pode ser substancial, particularmente para kernels de ordem superior. O \(\mbox{MISE}\) exato pode ser obtido por integração numérica das equações de viés e variância. Para dados normais, o \(\mbox{MISE}\) exato foi obtido para vários tamanhos de amostra com os kernels da Tabela 6.1, além do histograma.

    Os resultados são apresentados na Figura 6.4. As curvas individuais do \(\mbox{MISE}\) são mostradas graficamente em função de \(h/h^*\), de modo que o mínimo esteja centrado em 1. Essas figuras sugerem que, na maioria das situações práticas, kernels de ordem 2 e 4 são suficientes. O maior ganho em \(\mbox{MISE}\) é obtido ao passar do histograma para o kernel de ordem 2, com retornos decrescentes a partir daí. Kernels de ordem superior também parecem sensíveis ao excesso de suavização. Os kernels de ordem 8 são inferiores ao histograma se \(h>2h^*\).

    Figura 6.4: \(\mbox{MISE}\) exato usando kernels de ordem superior com dados normais para vários tamanhos de amostra. O \(\mbox{MISE}\) do histograma está incluído para referência.

    Esses kernels de ordem superior possuem componentes negativos. Eles serão chamados de “kernels negativos”, embora a expressão mais precisa seja “não não-negativos”. A introdução de kernels negativos proporciona uma melhoria no \(\mbox{MISE}\), mas ao custo de explicações especiais. Essa negatividade é particularmente incômoda em múltiplas dimensões, onde as regiões de estimação negativa podem estar espalhadas por todo o domínio.

    Estatísticos podem se sentir confortáveis ignorando tais características, mas deve-se ter cuidado em seu uso prático. Na prática, porções negativas da estimativa podem ser truncadas em 0. O truncamento introduz descontinuidades na derivada da estimativa, e a estimativa de densidade modificada agora integra para um valor ligeiramente superior a 1.

    Na Figura 6.5, as estimativas \(\mbox{ASH}\) da equação (5.5), utilizando os kernels de ordem 2 e 4 da Tabela 6.1 para os pesos da equação (5.6), são aplicadas aos dados da superfície do aço (Bowyer 1980; Silverman 1986). Os dados são medições, a partir de uma origem arbitrária, da altura real de uma superfície plana usinada em uma grade de 15.000 pontos.

    Figura 6.5: Estimativas de kernel positivas e negativas dos dados da superfície do aço. Os kernels utilizados foram \(K_2\) e \(K_4\) da Tabela 6.1.

    As larguras de banda foram selecionadas de forma que os valores na moda coincidissem. Este exemplo indica claramente a razão pela qual muitos estatísticos estão dispostos a usar kernels negativos com grandes conjuntos de dados.

    O problema da negatividade é bastante marginal, valor mínimo -0.00008, e há alguma melhoria na aparência da estimativa no pico, menor viés e variância simultaneamente.

    Um segundo problema introduzido pelo uso de kernels negativos envolve a escolha subjetiva da largura de banda. Com kernels negativos, as estimativas podem parecer grosseiras para valores moderadamente suavizados de \(h\), e não apenas para valores pequenos e pouco suavizados, como ocorre com kernels positivos.

    Essa dupla irregularidade pode ser um problema para iniciantes, especialmente considerando a promessa de estimativas “melhores” de ordem superior. Esse fenômeno é demonstrado na Figura 6.6 com os dados de queda de neve.

    Figura 6.6: Kernel \(K_4\) aplicado aos dados de queda de neve em Buffalo com três parâmetros de suavização. A estimativa de \(\mbox{ASH}\) é representada em forma de histograma.

    É interessante especular sobre os méritos relativos de usar um kernel de ordem superior em vez de um kernel de ordem inferior com um procedimento adaptativo. Um kernel adaptativo se comporta de maneira muito semelhante a um polígono de frequência adaptativo, como no Teorema 4.2. Há razões para acreditar que ambos têm sua utilidade.

    Assintoticamente, kernels de ordem superior superam os procedimentos adaptativos. A sensibilidade a erros na largura de banda ótima aumenta com a ordem do kernel, sugerindo que a adaptabilidade é mais importante com kernels de ordem superior.

    Considere, por exemplo, o \(\mbox{MSE}(h,x)\) exato para a densidade \(N(0,1)\) com um kernel normal. Uma expressão analítica para \(\mbox{MSE}(h,x)\) foi dada por Fryer (1976). Um gráfico de \(\mbox{MSE}\) versus \(h\) para vários valores de \(x\) entre 0 e 3 é mostrado na Figura 6.7. A melhor largura de banda global para \(n = 1000\) é \(h^* = 0.266\). No gráfico da esquerda, a melhor largura de banda, conforme \(x\) varia, está próxima desse valor, exceto quando \(x\approx 1\). Lembre-se de que o viés é uma função de \(f''(x)\) e a segunda derivada para uma densidade normal se anula em \(x =\pm 1\); portanto, a largura de banda local ótima será maior nesse ponto para uma variância reduzida.

    Ao mostrar graficamente em uma faixa mais ampla de larguras de banda, revela-se uma surpresa no gráfico da direita da Figura 6.7. Para os pontos \(x\) na cauda da distribuição, existem algumas larguras de banda muito grandes que resultam em \(\mbox{MSE}\)s surpreendentemente pequenos. A razão é que, ao calcular a média em uma vizinhança muito grande, o viés pode ser eliminado, ajustando-se artificialmente à densidade verdadeira.

    Figura 6.7: O erro quadrático médio exato em função da largura de banda \(h\) para \(f\sim N(0,1)\) e \(n = 1000\) para valores selecionados de \(x\) entre 0 e 3. A largura de banda globalmente ótima \(h = 0.266\) é indicada pela linha pontilhada vertical. O intervalo de largura de banda é expandido no quadro à direita.

    A utilidade potencial dessas grandes larguras de banda é discutida abaixo na Seção 6.6.4.1. Existe também a melhor largura de banda adaptativa “local” em torno de \(h = 0.266\), que é uma meta razoável. Schucany (1989) relatou algum sucesso nessa tarefa. Mas, no geral, uma boa estimativa adaptativa continua sendo uma tarefa difícil. A função \(\mbox{MSE}\) é notavelmente complexa, dada a simplicidade da densidade \(N(0,1)\).

    Vamos retornar à discussão sobre a abordagem de kernel de ordem superior. Embora tomar o limite dos kernels quando \(p\to\infty\) possa não parecer prudente, Davis (1975) investigou as propriedades de um kernel específico para “\(p=\infty\)”, a função “sinc” \(\sin(x)/(\pi x)\). Foi demonstrado que o \[ \mbox{MISE}^* = O\big(n/\log(n)\big)\cdot \]

    Marron and Wand (1992) examinaram o \(\mbox{MISE}\) de uma variedade de densidades e kernels mais complexos na família de misturas normais e calcularam os tamanhos de amostra necessários para justificar o uso de um kernel de ordem superior. Finalmente, kernels de ordem superior podem ser usados para estimar a derivada de uma densidade. No entanto, dada a aparência não monotônica de tais kernels, as estimativas da derivada provavelmente apresentarão artefatos de kernel e devem ser reservadas para situações com grande quantidade de dados.

    G. R. Terrell and Scott (1980), utilizando uma técnica de jackknife generalizada semelhante à de Schucany and Sommers (1977), propuseram um método alternativo para reduzir o viés. O método jackknife reduz o viés comparando dois estimadores. Na estimação de densidade, o procedimento envolve a construção da razão entre dois estimadores de kernel positivos com diferentes larguras de banda; por exemplo, \[ \tag{6.24} \widehat{f}(x)=\widehat{f}_h(x)^{4/3} \Big/ \widehat{f}_{2h}(x)^{1/3}\cdot \] O resultado decorre da aplicação do operador jackknife ao termo de viés \(O(h^2)\) em uma expansão em série de Taylor do logaritmo do viés.

    A esperança \(\mbox{E}_h\) da estimativa usual do kernel é \[ \mbox{E}_h=\mbox{E}\big(\widehat{f}_h(x) \big)=f(x)\left(1+\dfrac{c_2 h^2}{f(x)}+\dfrac{c_4 h^4}{f(x)}+\cdots \right), \] onde \(c_2=h^2 \sigma_K^2 f''(x)/2\), etc.

    Então, com a expansão de Taylor, \[ \begin{array}{rcl} \log \big( \mbox{E}_h\big) & = & \displaystyle \log\big(f(x)\big) + \dfrac{c_2 h^2}{f(x)}+ \dfrac{\big(c_4 f(x)-c_2^2/2 \big)h^4}{f(x)^2}+\cdots \\[0.8em] \log \big( \mbox{E}_{2h}\big) & = & \displaystyle \log\big(f(x)\big) + \dfrac{4 c_2 h^2}{f(x)}+\dfrac{\big( c_4 f(x)-c_2^2/2 \big) h^4}{f(x)^2}+\cdots \end{array} \]

    Então, \[ \dfrac{4}{3}\log\big( \mbox{E}_h\big) -\dfrac{1}{3}\mbox{E}_{2h} = \log\big(f(x)\big)+O(h^4), \] que, após a aplicação de exponenciais, sugere a forma dada na equação (6.24), veja o Exercício 14 para detalhes.

    Os autores mostram que \(h^* = 1.42\, \sigma \, n^{-1/9}\) no caso \(f = N(0,1)\) e \(K = N(0,\sigma^2)\). A estimativa resultante é não negativa e contínua, mas sua integral geralmente é ligeiramente maior que 1.

    Geralmente, exceder a taxa \(n^{-4/5}\) requer a violação de um dos dois pressupostos de uma função de densidade: não negatividade e massa de probabilidade total igual a 1.


    6.2.3.2 Kernels ótimos


    A estimativa de densidade do kernel herda todas as propriedades do seu kernel. Portanto, é importante notar que o kernel de Rosenblatt ingênuo é descontínuo, o kernel triangular \(\mbox{ASH}\) tem uma derivada descontínua e o kernel de Cauchy não possui momentos.

    Uma recomendação conservadora é escolher um kernel suave, claramente unimodal e simétrico em relação à origem. No entanto, formatos de kernel estranhos raramente são visíveis na estimativa final, exceto talvez nas caudas, devido a todas as médias.

    A questão de encontrar um kernel ótimo para estimativas não negativas foi considerada por Epanechnikov (1969); o mesmo problema variacional foi considerado por Bartlett (1963) e, em outro contexto, por Hodges and Lehmann (1956).

    A partir da equação (6.18), a contribuição do kernel para o \(\mbox{AMISE}\) ótimo é o seguinte fator adimensional: \[ \tag{6.25} \mbox{AMISE}^* \propto \big(\sigma_K \, R(K) \big)^{4/5} \cdot \]

    O problema de encontrar a densidade mais suave para o problema de largura de banda excessivamente suavizada é semelhante ao problema de minimizar (6.25), que pode ser escrito como \[ \min_K \; R(K) \qquad \mbox{restrito à} \qquad \sigma_K^2=\sigma\cdot \]

    A solução é uma versão escalonada do chamado kernel de Epanechnikov: \[ K_2^*(t) =\dfrac{3}{4}(1-t^2) \pmb{I}_{[-1,1]}(t)\cdot \] É interessante notar que o kernel ótimo possui suporte finito. O kernel ótimo não é diferenciável em \(t = \pm 1\).

    A variância de \(K_2^*\) é 1/5 e \(R(K_2^*) = 3/5\) em (6.25). Como o \(\mbox{AMISE}\) também é proporcional a \(n^{-4/5}\), outros kernels requerem (ver Exercício 8) \[ \tag{6.26} \dfrac{\sigma_K \, R(K)}{\sigma_{K_2^*} \, R(K_2^*)} = \dfrac{\sigma_K \, R(K)}{3/(5\sqrt{5})} \] vezes mais dados para atingir o mesmo AMISE que o kernel Epanechnikov ótimo.

    A Tabela 6.2 lista muitos kernels comumente usados e calcula sua eficiência relativa assintótica. O kernel ótimo mostra apenas uma melhoria modesta. Portanto, o kernel pode ser escolhido por outros motivos, facilidade de computação, diferenciabilidade, etc. sem grande preocupação com a perda de eficiência. É um tanto surpreendente que o popular kernel normal seja tão dispendioso. Dado o custo computacional do cálculo de exponenciais, é difícil recomendar o uso efetivo do kernel normal, exceto como ponto de referência.

    Tabela 6.2: Alguns kernels comuns e alguns incomuns e suas respectivas eficiências. Todos os kernels tem suporte em \([-1,1]\), a menos que indicado de outra forma.

    Na última parte da tabela, alguns kernels absurdos são listados para ilustrar que uma ineficiência muito grande ainda é menor que 2. Pelo Teorema 4.1, o estimador do polígono de frequência (FP) pertence ao mesmo grupo dos estimadores de kernel positivo.

    As entradas na tabela referentes ao FP foram obtidas pela correspondência das expressões \(\mbox{AMISE}\) nos Teoremas 4.1 e 6.1. A conclusão é que o FP está de fato na mesma classe, mas é ineficiente. Finalmente, observe que o kernel limite do histograma deslocado médio (triângulo isósceles) é superior ao kernel limite do polígono de frequência deslocado médio (FP indiferente), embora o próprio FP seja superior ao histograma.

    As funções de densidade Beta simétricas, quando transformadas para o intervalo \((-1,1)\) de modo que a média seja 0, são uma escolha útil para uma classe de kernels: \[ \tag{6.27} K_k(t)=\dfrac{(2k+1)!!}{2^{k+1}k!} (1-t^2)_+^k, \] onde a notação de fatorial duplo significa \((2k+1)!! = (2k+1)(2k-1)\cdots 5\cdot 3\cdot 1\). Os kernels de Epanechnikov e biweight pertencem a esta classe. O mesmo ocorre com a densidade normal quando \(k\to\infty\) (ver Exercício 18).

    A busca por kernels ótimos de alta ordem é bastante diferente e não tão frutífera. Suponha que \(K_4^*\), \(K_6^*\), … sejam os kernels ótimos de ordem 4, ordem 6, …, respectivamente. Pelo Teorema 6.2, a contribuição do kernel para o \(\mbox{AMISE}^*\) para kernels de ordem 4 é a seguinte quantidade não negativa e adimensional: \[ \tag{6.28} \mbox{AMISE}^* \propto \big( R(K)^8 \mu_4^2\big)^{1/9}\cdot \] Considere o seguinte kernel de quarta ordem, que é uma mistura de \(K_4^*\) e \(K_6^*\): \[ K_\epsilon(t)=\epsilon K_4^*(t)+(1-\epsilon) K_6^*(t), \qquad 0\leq \epsilon\leq 1\cdot \]

    \(K_\epsilon\) possui rugosidade finita, mas seu quarto momento se anula quando \(\epsilon\to 0\). Portanto, o quarto momento de \(K_4^*\) deve ser 0 e o critério em (6.28) é igual a 0 na solução. Mas, por definição, \(K_4^*\) deixa de ser um kernel de quarta ordem. Como muitos kernels têm quarto “momento” igual a zero, mas rugosidade finita, o limite inferior de 0 é atingido por muitos kernels, nenhum dos quais é de ordem 4 ou interessante no sentido de Epanechnikov.

    Obviamente, o \(\mbox{AMISE}\) não seria de fato 0, mas envolveria termos de ordem superior. A escolha entre kernels de ordem superior é bastante complexa e é difícil estabelecer diretrizes. Na prática, os métodos de segunda e quarta ordem são provavelmente os que mais se devem considerar, já que kernels além da ordem 4 proporcionam pouca redução adicional no \(\mbox{MISE}\).

    Para amostras muito grandes, onde os métodos de ordem superior proporcionam uma redução fracionária substancial, o \(\mbox{MISE}\) absoluto pode já ser tão pequeno que qualquer vantagem prática se perde.

    A recomendação geral para a escolha de um kernel com base nessas observações é escolher um kernel simétrico que seja um polinômio de baixa ordem. Gasser, Müller, and Mammitzsch (1985) desenvolveram uma hierarquia suave de kernels de ordem superior. Eles mostram que seus kernels são ótimos, mas em um sentido diferente: esses kernels têm rugosidade mínima sujeita a um número fixo de mudanças de sinal no kernel. Tal justificativa não justifica realmente o rótulo “ótimo” no sentido de Epanechnikov. No entanto, eles forneceram uma fórmula valiosa para kernels polinomiais de baixa ordem apropriados para várias combinações da ordem \(p\) do kernel e para a \(r\)-ésima derivada da densidade, \[ K_{(k,p,r)}^*=\sum_{i=0}^{k+2(r-1)} \lambda_i \, t^i \, \pmb{I}_{[-1,1]}(t), \] onde \(k\geq p+2\) e \((k,p)\) são ambos ímpar or ambos pares, com \[ \lambda_i = \dfrac{(-1)^{(i+p)/2} (k+p+2r)! \, (k-p)(k+2r-i) (k+i)!}{i! \, (i+p+1) 2^{2(k+r)+1} \left(\frac{k-p}{2} \right)! \left(\frac{k+p+2r}{2} \right)!\left(\frac{k+2r-i}{2} \right)!\left(\frac{k+i}{2} \right)!}, \] se \(k+i\) for ímpar, e 0 caso contrário (ver Müller 1988).

    Dada a disponibilidade de programas de manipulação simbólica, provavelmente é suficiente resolver o conjunto de equações lineares que regem a aplicação específica de forma simples.

    Uma abordagem de “kernel projetado” permite a adição de quaisquer condições lineares e resulta em um novo kernel de grau polinomial mais elevado. Novamente, a escolha do kernel não é um fator crítico.


    6.2.3.3 Kenels equivalentes


    Por uma série de razões, não existe um único núcleo que possa ser recomendado para todas as circunstâncias. Um candidato viável é o núcleo normal; no entanto, ele é relativamente ineficiente e possui suporte infinito.

    O kernel de Epanechnikov ótimo não é continuamente diferenciável e não pode ser usado para estimar derivadas. Na prática, a capacidade de alternar entre diferentes kernels sem ter que reconsiderar o problema de calibração a cada mudança é conveniente. Essa tarefa é fácil de realizar, mas apenas para kernels da mesma ordem.

    Como Scott (1976) observou, se \(h_1\) e \(h_2\) são parâmetros de suavização a serem usados com os kernels \(K_1\) e \(K_2\), respectivamente, então o Teorema 6.1 implica que assintoticamente \[ \tag{6.29} \dfrac{h_1^*}{h_2^*}=\left(\dfrac{R(K_1)/\sigma_{K_1}^4}{R(K_2)/\sigma_{K_2}^4} \right)^{1/5} = \dfrac{\sigma_{K_1}}{\sigma_{K_2}}\left( \dfrac{\sigma_{K_1} R(K_1)}{\sigma_{K_2} R(K_2)}\right)^{1/5}\cdot \] A Tabela 6.3 apresenta um resumo dos fatores para larguras de banda de suavização equivalentes entre kernels populares.

    Tabela 6.3: Fatores para suavização equivalente entre kernels populares. Para ir de \(h_1\) para \(h_2\), multiplique \(h_1\) pelo fator na tabela, na linha rotulada como \(K_1\) e na coluna rotulada como \(K_2\).

    O termo entre parênteses à direita na equação (6.29) é a razão entre quantidades adimensionais para cada kernel. Essas quantidades \(\sigma_K R(K)\) são quase iguais entre si, como pode ser observado na equação (6.26) e na Tabela 6.2.

    Assim, o termo entre parênteses na equação (6.29) pode ser definido como 1, de modo que a tarefa de escolher parâmetros de suavização equivalentes para diferentes kernels possa ser realizada por meio de um escalonamento de acordo com os desvios padrão: \[ \tag{6.30} \mbox{Redimensionamento equivalente do kernel: } \quad h_2^*\approx \dfrac{\sigma_{K_1}}{\sigma_{K_2}} h_1^*\cdot \]

    Por exemplo, o fator “exato” para passar de uma largura de banda normal para uma largura de banda triweight na Tabela 6.3 é 2.978, enquanto a regra aproximada (6.30) fornece o fator 3.

    O escalonamento de largura de banda equivalente fornece estimativas quase idênticas não apenas para parâmetros de suavização ótimos, mas também para valores não ótimos. Esse reescalonamento é frequentemente usado ao calcular e mostrar graficamente uma estimativa de kernel biweight, mas usando um parâmetro de suavização derivado de uma regra de validação cruzada de kernel normal.

    Se todos os kernels fossem apresentados com variâncias iguais, nenhuma alteração nos parâmetros de suavização seria necessária. No entanto, seria extremamente difícil lembrar as fórmulas desses kernels, pois as variâncias estariam incorporadas às formas dos kernels. Em suma, parece mais fácil escrever kernels em uma forma parcimoniosa. Marron and Nolan (1988) propuseram escalonar todos os kernels para sua forma “canônica”, que possui larguras de banda equivalentes a um kernel normal. Esta proposta difere ligeiramente da proposta de reescalonamento da variância, uma vez que a rugosidade do kernel também é levada em consideração.

    Para kernels de ordem superior \(K_{p_{_1}}\) e \(K_{p_{_2}}\) de ordem \(p\), um reescalonamento semelhante decorre do Teorema 6.3 com base no “momento” de ordem superior apropriado (ver Exercício 20): \[ \tag{6.31} \dfrac{h_{p_{_1}}^*}{h_{p_{_2}}^*}=\left(\dfrac{\mu_{p_{_2}}}{\mu_{p_{_1}}} \right)^{1/p} \left(\left(\dfrac{\mu_{p_{_1}}}{\mu_{p_{_2}}} \right)^{1/p} \dfrac{R(K_{p_{_1}})}{R(K_{p_{_2}})} \right)^{1/(2p+1)}\approx \left(\dfrac{\mu_{p_{_2}}}{\mu_{p_{_1}}} \right)^{1/p}, \] visto que a quantidade entre parênteses é adimensional e aproximadamente igual a 1 para a maioria dos kernels simétricos.

    Quando \(p = 2\), as equações (6.31) e (6.30) coincidem. Alguns dos momentos de ordem superior são negativos, mas sua razão é positiva em (6.31).


    6.2.3.4 Kernels de ordem superior e projeto de kernel


    Entre kernels de diferentes ordens, não existe uma noção similar para a escolha de larguras de banda que proporcionem “suavização equivalente”. Além disso, a ausência de um verdadeiro “kernel ótimo” além da ordem 2 é preocupante. Nesta seção, kernels de ordem superior são reintroduzidos a partir de dois outros pontos de vista.

    Os resultados têm algumas implicações curiosas para a seleção de largura de banda e o projeto de kernel, e sugerem maneiras pelas quais trabalhos futuros podem ser úteis.

    A primeira abordagem alternativa recorre ao argumento da análise numérica introduzido inicialmente por Rosenblatt (1956). Utilizando aproximações de diferenças finitas progressivas e centrais para a derivada, foi demonstrado na Seção 6.1.1 que os núcleos equivalentes são os núcleos de ordem 1 e ordem 2 \(U(0,1)\) e \(U(-0.5, 0.5)\), respectivamente.

    Utilizando a notação \[ \Delta F_n(x,c) = F_n(x+c)-F_n(x-c), \] o estimador de Rosenblatt de dois pontos com kernel \(U(-1,1)\) é \(\widehat{f}_2(x) = \Delta F_n (x,h)/(2h)\).

    Considere os seguintes estimadores de derivadas de 4, 6 e 8 pontos (veja o Exercício 21): \[ \tag{6.32} \begin{array}{rcl} \widehat{f}_4(x) & = & \displaystyle \dfrac{8\Delta F_n(x,h)-\Delta F_n(x,2h)}{12h},\\[0.8em] \widehat{f}_6(x) & = & \displaystyle \dfrac{45 \Delta F_n(x,h)-9\Delta F_n(x,2h)+ \Delta F_n(x,3h)}{60h},\\[0.8em] \widehat{f}_8(x) & = & \displaystyle \dfrac{224\Delta F_n(x,h)-56 \Delta F_n(x,2h)+\frac{32}{3}\Delta F_n(x,3h)-\Delta F_n(x,4h) }{280h}\cdot \end{array} \] Os vieses para esses três estimadores são \(-h^4 f^{(4)}(x)/30\), \(-h^6 f^{(6)} (x)/140\) e \(-h^8 f^{(8)(x)/630\), respectivamente.

    Os kernels equivalentes podem ser visualizados mostrando \(\widehat{f}_r(x)\) com um ponto de dados em 0 com \(h = 1\), como na Figura 6.8.

    Figura 6.8: Núcleos (kernels) boxcar ou ingênuos de ordem superior baseados em diferenças finitas.

    Na Figura 6.9, a estimativa do kernel ingênuo \(p = 2\) é mostrada para 320 níveis de colesterol de pacientes com doença cardíaca. Como \(\widehat{\sigma}= 43\), a largura de banda escolhida foi \(h = 25\), que é a regra de largura de banda equivalente para um kernel \(U(-1,1)\) à regra de referência normal (6.19).

    Figura 6.9: Kernels (núcleos) ingênuos de ordem superior aplicados aos dados de colesterol de homens doentes \((n = 320)\). As larguras de banda são indicadas pelas linhas horizontais.

    Especificamente, a largura de banda equivalente é calculada como \[ h\approx 3^{1/2} \big(1.06\times 43\times 320^{-1/5}\big), \] onde \(3^{1/2}\) é o desvio padrão do kernel \(U(-1,1)\).

    As larguras de banda para os kernels de ordem $p = $ quarta, sexta e oitava foram escolhidas de forma que os níveis nas modas fossem iguais. Os lóbulos laterais negativos são facilmente visíveis.

    As estimativas de segunda e quarta ordem são substancialmente diferentes. A linha horizontal mostra a largura de banda \(h\), que é metade da largura do lóbulo central do kernel. Assim, para \(K_8\), o suporte é oito vezes maior que \(h\); a influência de cada ponto de dados se estende por 216 unidades (mg/dl) em ambas as direções. Para os kernels de ordem 2, 4 e 6, a extensão é de 25, 92 e 153, respectivamente. Kernels de ordem superior, mesmo com suporte finito, não são locais.

    Uma observação empírica é que, quando um kernel de ordem superior não consegue reduzir ainda mais o \(\mbox{ISE}\), as estimativas de densidade de ordem superior ótimas são todas muito semelhantes, exceto por um pequeno ruído local devido à rugosidade adicional nas caudas do kernel.

    Essa semelhança ocorre quando o lóbulo central do kernel, no intervalo \((-h,h)\), permanece com largura fixa mesmo com o aumento da ordem do kernel. A sequência de larguras de banda é 25, 46, 51 e 54 para os dados de colesterol. Assim, \(p = 4\) parece uma escolha plausível para esses dados, possivelmente com uma largura de banda menor. À medida que \(p\) aumenta, os lóbulos negativos têm uma tendência indesejável de crescer e se espalhar.

    Alternativamente, as larguras de banda para os kernels de ordem superior ingênuos poderiam ter sido escolhidas para aumentar a estimativa na moda por uma estimativa do viés ali; por exemplo, para o kernel de ordem 2, \[ \tag{6.33} \mbox{Viés}\big(\widehat{f}(x) \big)=\dfrac{1}{2}h^2 \sigma_K^2 f''(x)\cdot \] Se a estimativa resultante, com o kernel de ordem 4, for muito grosseira ou se não for possível encontrar uma largura de banda que proporcione o aumento desejado, então uma conclusão razoável é que a ordem máxima viável do kernel foi excedida. A largura de banda resultante para o kernel de ordem superior ingênuo pode ser transformada por (6.31) em um kernel de ordem superior mais suave na Tabela 6.1.

    A segunda introdução alternativa a kernels negativos sugere que o problema da seleção da largura de banda é ainda mais complexo do que aparenta. Partindo de uma estimativa de kernel positiva e da estimativa de viés fornecida anteriormente, a ideia é estimar e remover o viés pontualmente usando uma estimativa de kernel da segunda derivada.

    Para maior clareza, uma largura de banda \(g\) possivelmente diferente é usada para estimar \(f''(x)\): \[ \tag{6.34} \widehat{f}''(x)=\dfrac{1}{ng^3}\sum_{i=1}^n K'' \left( \dfrac{x-x_i}{g}\right)\cdot \] Considere a estimativa do kernel “corrigida para o viés”, que é obtida combinando as equações (6.33) e (6.34): \[ \tag{6.35} \widehat{f}(x)=\dfrac{1}{nh}\sum_{i=1}^n K\left(\dfrac{x-x_i}{h} \right)-\dfrac{1}{2}h^2 \sigma_k^2 \times \dfrac{1}{ng^3} \sum_{i=1}^n K'' \left( \dfrac{x-x_i}{g}\right)\cdot \] Como \(\widehat{f}''(x)\) integra para 0, essa estimativa de kernel modificada ainda integra para 1. Por inspeção, (6.35) está na forma de um estimador de kernel com kernel \[ \tag{6.36} K_{h,g}(t)=K_h(t)-\dfrac{h^2 \sigma_K^2}{2g^2} K_g''(t)\cdot \] Um cálculo direto verifica que este novo kernel tem um segundo momento nulo, de modo que este procedimento construtivo resulta, na verdade, em um kernel de ordem 4 (veja o Exercício 22).

    No entanto, esta abordagem sugere que \(g = h\), uma vez que a largura de banda para a segunda derivada deve ser maior do que para a própria estimativa de densidade. Se assim for, um kernel de ordem superior bem construído deve “expandir-se” lentamente à medida que o tamanho da amostra aumenta. Dadas as mudanças relativamente lentas nas larguras de banda, pode não haver nenhuma melhoria prática em relação a permitir \(g = h\).


    6.2.3.5 Núcleos de fronteira


    Quando a densidade desconhecida é descontínua, as estimativas por núcleo sofrem a mesma perda drástica de eficiência \(\mbox{MISE}\) que o polígono de frequência. No caso em que a localização da descontinuidade é conhecida, uma solução elegante está disponível.

    Sem perda de generalidade, suponha que a descontinuidade ocorra em zero e que a densidade se anule para \(x < 0\). Uma análise cuidadosa do argumento teórico que levou à eliminação do termo de viés \(O(h)\) para uma estimativa por núcleo revela que o requisito crítico é que o núcleo satisfaça \[ \int t K(t)\mbox{d}t = 0\cdot \] O requisito fundamental não é que o núcleo seja simétrico; a simetria é apenas uma maneira simples pela qual o kernel pode satisfazer a equação integral.

    A tarefa, portanto, é projetar kernels de suporte finito para uso com amostras na região de fronteira \(x_i\in (0,h)\). Para que a estimativa do kernel se anule para \(x < 0\), o kernel para um ponto de amostra \(x_i\in [0,h)\) deve cobrir o intervalo \([0,x_i+h)\) em vez do intervalo \((x_i-h, x_i+h)\).

    Como o intervalo \([0, x_i+h)\) é mais estreito que \(2h\), a rugosidade dos kernels e, portanto, a \(\mbox{IV}\), aumentará consideravelmente. Em um contexto de regressão, Gasser and Müller (1979) e Rice (1984) sugeriram o uso do intervalo mais amplo \((0, 2h)\) para cada \(x_i\in [0,h)\). Esta sugestão é equivalente a escolher kernels com suporte no intervalo \((c, c+2)\), para \(-1\leq c \leq 0\), que satisfaçam \[ \int_c^{c+2} t K(t)\mbox{d}t = 0\cdot \]

    Uma tentativa de permitir que a largura do intervalo varie para alcançar “suavização equivalente” usando a regra completa (6.29) parece fadada ao fracasso, pois um intervalo mais amplo geralmente não pode ser encontrado com suavização equivalente. Assim, a escolha simples de \(2h\) é um compromisso razoável.

    Descreve-se um kernel de fronteira projetado. Suponha que o kernel de fronteira desejado seja uma modificação do kernel biweight comum, com propriedades semelhantes no ponto final direito \(x = c+2\). Isso sugere a análise de um kernel de fronteira projetado da forma \[ K_c(t)=\big(c_1+c_2(t-c)^2 \big)\times \big(t-(c+2) \big)^2, \qquad -1\leq c\leq 0, \] onde as constantes \(c_1\) e \(c_2\) são determinadas pelas restrições \[ \int K_c (t)\mbox{d}t = 1 \qquad \mbox{e} \qquad \int t K_c(t) \mbox{d}t = 0\cdot \] A forma para \(K_c(\cdot)\) garante que as duas restrições sejam lineares nas incógnitas \(c_1\) e \(c_2\).

    A resolução dessas equações resulta em \[ \tag{6.37} K_c(t)=\dfrac{3}{4}\left( (c+1)-\dfrac{5}{4}(1+2c)(t-c)^2 \right)\times \big(t-(c+2) \big)^2 \pmb{I}_{[c,c+2]}(t)\cdot \]

    A Figura 6.10 mostra alguns exemplos desses kernels de duas maneiras diferentes. Primeiro, os kernels são mostrados centrados na amostra (considerada a origem). Segundo, os kernels são mostrados como seriam colocados sobre as amostras, de modo que comecem em zero. Observe que se \(x_i\in [0,h)\), então o kernel \(K_c\) com \(c = (0-x_i)/h\) deve ser usado em vez do kernel biweight comum.

    Figura 6.10: (Quadro esquerdo) Exemplos dos kernels de fronteira “flutuantes” \(K_c(t)\), onde \(-1 < c < 0\). (Quadro direito) Assumindo a fronteira \(x\geq 0\), cada kernel \(K_c(t)\) é desenhado centrado no ponto de dados, \(x_i = -c\), que é indicado pela linha vertical tracejada.

    Claramente, a Figura 6.10 indica que esses são kernels de fronteira flutuantes, o que significa que o valor do kernel flutua na fronteira esquerda. Um exemplo do uso desses kernels com uma amostra de 100 pontos da densidade exponencial negativa ilustra a eficácia dos kernels de fronteira flutuantes, veja a Figura 6.11.

    Figura 6.11: (Quadro esquerdo) Exemplo com dados exponenciais negativos, com e sem modificação de fronteira para \(n = 100\) e \(h = 0.93\). Os kernels de fronteira “flutuante” e “zero” são definidos nas equações (6.37) e (6.38), respectivamente. (Quadro direito) Exemplo com densidade \(Beta(3,9)\) em uma vizinhança de 0 para \(n = 100\) e \(h = 0.11\).

    Observe como a estimativa do kernel biweight não modificado é bastante enviesada no intervalo \((0, h)\) e se estende para a região \(x < 0\). No entanto, sem qualquer verificação ou indicação de um problema de fronteira, a estimativa do kernel não modificado parece bastante suave. Essa suavidade deve servir como um alerta para verificar erros resultantes da falta de conhecimento prévio da existência de um problema de fronteira.

    Para a densidade exponencial negativa, a teoria assintótica é válida se a rugosidade for calculada no itervalo \((0,\infty)\): \[ \int_0^\infty f''(x)^2 \mbox{d}x=1/2, \] de modo que \(h^*=(70/n)^{1/5}\) para o kernel biweight pelo Teorema 6.1.

    A estimativa de densidade dos dados exponenciais negativos, mostrada como uma linha tracejada grossa na Figura 6.11, foi construída usando o kernel de fronteira zero, \[ \tag{6.38} K_c^0(t) = \dfrac{15}{16}(t-c)^2 (2+c-t)^2 \big((7c^2+14c + 8)-7 t (c+1) \big)\pmb{I}_{[c,c+2]}(t)\cdot \]

    Este kernel biweight modificado foi projetado com a restrição adicional de que o kernel de fronteira e sua derivada se anulem na fronteira esquerda \(x = c\), em vez de flutuarem como antes. Essa modificação foi realizada na fase de projeto, incluindo o fator \((t-c)^2\), que garante que o kernel e sua derivada se anulem no ponto final esquerdo \(t = c\).

    A Figura 6.12 exibe esses kernels de duas maneiras. Claramente, o kernel é inadequado para dados exponenciais negativos, induzindo uma oscilação grande, porém suave. No entanto, o kernel de fronteira zero é apropriado para dados da densidade \(Beta(3,9)\); veja a Figura 6.11, que mostra a aplicação desses kernels a uma amostra de 100 pontos \(Beta(3,9)\) na vizinhança da origem. Para essa densidade, \(R(f'')\approx 24.835\) e \(h^* \approx (0.269/n)^{1/5}\). A estimativa comum se estende para a região negativa e a estimativa de kernel “flutuante” faz jus ao seu nome.

    Figura 6.12: Exemplos de kernels de fronteira “zero”, como na Figura 6.10.

    Embora os kernels de fronteira possam ser muito úteis, existem problemas potencialmente sérios com dados reais. Há um número infinito de kernels de fronteira que refletem o espectro de possíveis restrições de projeto, e esses kernels não são intercambiáveis.

    Artefatos graves podem ser introduzidos por qualquer um deles em situações inadequadas. É necessário um exame muito cuidadoso para evitar ser vítima do kernel de fronteira escolhido. Artefatos podem, infelizmente, ser introduzidos pela escolha do intervalo de suporte para o kernel de fronteira. Pouco se sabe sobre a melhor maneira de evitar essa situação, mas a solução de Rice-Müller parece ser a melhor dentre várias alternativas possíveis que foram tentadas.

    Finalmente, embora os kernels de fronteira pareçam cobrir uma ampla região da estimativa de densidade, o efeito geralmente se limita ao intervalo \((0, h/2)\) para estimativas suavizadas adequadamente.

    Alguns dados podem exigir um kernel de fronteira diferente em cada extremidade. Por exemplo, os 857 tempos mais rápidos da Maratona Tenneco de Houston, de 20 de janeiro de 1991, foram registrados como uma razão em relação ao tempo do vencedor. Claramente, espera-se que as características na densidade sejam diferentes nas duas fronteiras. Uma estimativa foi construída usando \(K_c^0\) à esquerda e o \(K_c\) flutuante à direita com \(h = 0.05\) (veja a Figura 6.13). O aglomerado entre os líderes é real; no entanto, a protuberância extra à direita parece ser mais um artefato.

    Figura 6.13: Estimativa de densidade dos 857 tempos mais rápidos na Maratona Tenneco de Houston de 1991. Os dados representam a razão entre o tempo do corredor e o tempo do líder da prova. Diferentes kernels de limite foram utilizados em cada extremo. Um histograma é apresentado para comparação.

    Reflexão acerca da técnica de fronteira: Existe uma técnica mais conservadora que pode substituir o kernel “flutuante”. Se os dados forem não negativos e a descontinuidade estiver em \(x = 0\), uma estimativa de kernel comum é calculada, mas nos dados aumentados \[ (-x_n \cdots,-x_1, x_1,\cdots, x_n)\cdot \]

    A estimativa final é obtida dobrando-se essa estimativa para \(x\geq 0\). A largura de banda deve ser baseada no tamanho da amostra \(n\) e não em \(2n\). Essa técnica evita as armadilhas dos kernels de fronteira negativos, mas geralmente apresenta consistência de ordem inferior, veja o Exercício 4 no contexto do polígono de frequência.


    6.3 Propriedades teóricas: caso multivariado


    A análise teórica dos estimadores de kernel multivariados é a mesma que a dos polígonos de frequência, exceto por alguns detalhes. A discussão inicial será limitada aos estimadores de densidade de kernel produto. A análise geral de kernel será considerada posteriormente.


    6.3.1 Kernels produto


    A forma geral de um estimador de kernel (núcleo) produto é dada por \[ \tag{6.39} \widehat{f}(\pmb{x})=\dfrac{1}{n\, h_1 \, \cdots \,d_d} \sum_{i=1}^n \left( \prod_{j=1}^d K\left(\dfrac{x_i-x_j}{h_j} \right)\right)\cdot \] O mesmo kernel (univariado) é usado em cada dimensão, mas com um parâmetro de suavização (possivelmente) diferente para cada dimensão. Os dados \(x_{ij}\) provêm de uma matriz \(n\times d\).

    A estimativa é definida pontualmente, onde \(\pmb{x} = (x_1,\cdots, x_d)^\top\). Geometricamente, a estimativa coloca uma massa de probabilidade de tamanho \(1/n\) centrada em cada ponto amostral, exatamente como no caso univariado. Lembre-se de que a forma limite do \(\mbox{ASH}\) multivariado ingênuo é um estimador de kernel triangular de produto. Vários kernels produto bivariados são exibidos na Figura 6.14.

    Figura 6.14: Exemplos de kernel produto para quatro kernels.

    Considere o viés pontual do estimador multivariado. Claramente, \[ \tag{6.40} \begin{array}{rcl} \mbox{E}\big( \widehat{f}(\pmb{x})\big) & = & \displaystyle \mbox{E}\left( \prod_{j=1}^d K\left(\dfrac{x_i-x_j}{h_j} \right)\right) = \int_{\mathbb{R}^d} \prod_{j=1}^d K\left(\dfrac{x_i-t_j}{h_j} \right) f(\pmb{t})\mbox{d}\pmb{t}\\[0.8em] & = & \displaystyle \int_{\mathbb{R}^d} \prod_{j=1}^d K\left(\omega_j\right) f(x_1-h_1\omega_1,\cdots,x_n-h_n\omega_n)\mbox{d}\pmb{\omega}\\[0.8em] & \approx & \displaystyle \int_{\mathbb{R}^d} \prod_{j=1}^d K\left(\omega_j\right)\left( f(\pmb{x})-\sum_{r=1}^d h_r\omega_r f_r(\pmb{\omega}) +\sum_{r,s=1}^d \dfrac{h_r h_s}{2} \omega_r\omega_s f_{rs}(\pmb{x})\right) \mbox{d}\pmb{\omega}\\[0.8em] & = & \displaystyle f(\pmb{x})+\dfrac{1}{2}\sigma_K^2 \sum_{j=1}^d h_j^2 f_{ij}(\pmb{x})+O(h^4)\cdot \end{array} \]

    Como antes, os termos de viés \(O(h)\) desaparecem se os kernels univariados tiverem média zero. Da mesma forma, os termos \(h_r h_s\) desaparecem (veja o Exercício 24). Segue-se que a integral do viés quadrático é como dado no Teorema 6.4. A variância pontual é \[ \dfrac{f(\pmb{x})R(K)^d}{n \, h_1\, h_2\cdots \, h_n}, \] da qual a integral da variância segue.


    Teorema 6.4:

    Para um estimador produto de kernel multivariado, os componentes do \(\mbox{AMISE}\) são \[ \tag{6.41} \begin{array}{rcl} \mbox{AISB} & = & \displaystyle\dfrac{1}{4}\sigma_K^4 \left(\sum_{i=1}^d h_i^4 R(f_{ij})+\sum_{i\neq j=1}^d h_i^2 \, h_j^2 \int_{\mathbb{R}^d} f_{ii}(\pmb{x}) f_{jj}(\pmb{x})\mbox{d}\pmb{x} \right), \\[0.8em] \mbox{AIV} & = & \displaystyle \dfrac{R(K)^d}{n\, h_1 \, h_2 \cdots \, h_d} -\dfrac{R(f)}{n}+O(\pmb{h}/n)\cdot \end{array} \]


    Demonstração. Scott (2015).


    A ordem dos parâmetros de suavização ótimos é exatamente a mesma que para o FP multivariado: \(h_i^∗ = O(n^{-1/(4+d)})\) e \(\mbox{AMISE}^∗ = O(n^{-4/(4+d)})\).

    É um pequeno inconveniente, mas não existe uma expressão analítica geral para os parâmetros de suavização ótimos, exceto como solução para \(d\) equações não lineares. Soluções podem ser encontradas em casos especiais, por exemplo, quando \(d\leq 2\) ou se \(h_i = h\) para todo \(i\).

    Por exemplo, com dados bivariados normais gerais e um kernel normal, uma integração direta mostra que \[ \begin{array}{rcl} R(f_{11}) & = & \displaystyle \dfrac{3}{16\pi (1-\rho^2)^{5/2} \, \sigma_1^2\sigma_2}, \\[0.8em] R(f_{12}) & = & \displaystyle \dfrac{3}{16\pi (1-\rho^2)^{5/2} \, \sigma_1 \sigma_2^5}, \\[0.8em] \displaystyle \int_{\mathbb{R}^2} f_{11}(x_1) f_{22}(x_2) \mbox{d}x_1 \mbox{d}x_2 & = & \dfrac{1+2\rho^2}{16\pi (1-\rho^2)^{5/2} \, \sigma_1^3 \sigma_2^3}\cdot \end{array} \]

    Pelo Teorema 6.4, o \(\mbox{AMISE}\) é minimizado quando \[ \tag{6.42} \begin{array}{rcl} h_i^* & = & \displaystyle \sigma_i (1-\rho^2)^{5/12} (1+\rho^2/2)^{-1/6} \, n^{-1/6} \\[0.8em] & \approx & \displaystyle \sigma_i (1-\rho^2/2 -\rho^4/16 - \cdots)\, n^{-1/6}, \qquad i=1,2, \end{array} \] para os quais \[ \mbox{AMISE}^* = \dfrac{3}{8\pi} \dfrac{1}{\sigma_1\sigma_2} \dfrac{1}{(1-\rho^2)^{5/6}}(1+\rho^2/2)^{1/3}\, n^{-2/3}\cdot \]

    Observe que o \(\mbox{AMISE}\) diverge para o infinito quando os dados são perfeitamente correlacionados, a verdadeira maldição da dimensionalidade. Em comparação com outras estimativas bivariadas, se \(\rho = 0\) e \(\sigma_i = 0\), então o \(\mbox{AMISE}\) bivariado é igual a 1/400 quando \(n = 302\). Um FP bivariado e um histograma requerem \(n = 557\) e \(n = 4244\), respectivamente.

    Um segundo exemplo de um caso especial é a distribuição normal multivariada, onde todas as variáveis são independentes. Se um kernel normal for usado, um breve cálculo com o Teorema 6.4 fornece a \[ \tag{6.43} \mbox{Regra de referência normal: } \; h_i^* =\left( \dfrac{4}{d+2}\right)^{1/(d+4)} \sigma_i n^{-1/(d+4)}\cdot \]

    À medida que a dimensão \(d\) varia, a constante na equação (6.43) varia no intervalo (0.924, 1.059), com um limite igual a 1. A constante é exatamente 1 no caso bivariado e mínima quando \(d = 11\).

    Portanto, uma regra amostral e fácil de lembrar é \[ \tag{6.44} \mbox{Regra de Scott em } \mathbb{R}^d: \; \widehat{h}_i=\widehat{\sigma}_i n^{-1/(d+4)}\cdot \]

    Para outros kernels, o parâmetro de suavização equivalente pode ser obtido dividindo-se pelo desvio padrão desse kernel. Assim como em uma dimensão, essas fórmulas podem ser usadas no lugar de valores de sobresuavização mais precisos, visto que dados normais independentes são muito suaves. Qualquer estrutura especial exigirá larguras de banda mais estreitas.

    Por exemplo, a modificação baseada na assimetria e curtose em 1 é idêntica aos fatores para o polígono de frequência na Seção 4.1.2. Se os dados não forem de posto completo, os métodos de kernel terão um desempenho ruim. Técnicas de redução de dimensionalidade serão consideradas no Capítulo 7.


    6.3.2 \(\mbox{MISE}\) de kernel multivariado geral


    Na prática, recomenda-se o uso de kernels produto. No entanto, para diversos estudos teóricos, serão necessários kernels multivariados gerais. Esta seção apresenta um breve resumo desses estudos.

    O estimador de kernel multivariado geral incluirá não apenas uma densidade multivariada arbitrária como kernel, mas também uma transformação linear arbitrária dos dados.

    Seja \(H\) uma matriz não singular \(d\times d\) e \(K \, : \, \mathbb{R}^d\to \mathbb{R}^1\) um kernel que satisfaça as condições abaixo. Então, o estimador de kernel multivariado geral é \[ \tag{6.45} \widehat{f}(\pmb{x})=\dfrac{1}{n|H|} \sum_{i=1}^n K\big( H^{-1}(\pmb{x}-\pmb{x_i})\big)\cdot \]

    Deve ficar evidente a partir da equação (6.45) que a transformação linear \(H\) pode ser incorporada à definição do kernel. Por exemplo, é equivalente escolher \(K\) como \(N(\pmb{0}_d, \Sigma)\) com \(H = \pmb{I}_d\), ou escolher \(K\) como \(N(\pmb{0}_d, \pmb{I}_d)\) com \(H = \Sigma^{1/2}\) (ver Exercício 25).

    Assim, é possível escolher um kernel multivariado com uma estrutura de covariância simples sem perda de generalidade. Não será suficiente, contudo, considerar apenas kernels produto, pois isso limitaria a discussão a kernels multivariados independentes, e não apenas não correlacionados, e a kernels com suporte em uma região retangular.

    A partir de agora, assumiremos que o kernel multivariado satisfaz três condições de momento, observe que estas são equações matriciais: \[ \tag{6.46} \begin{array}{rcl} \displaystyle \int_{\mathbb{R}^d} K(\pmb{\omega})\mbox{d}\pmb{\omega} & = & \displaystyle 1, \\[0.8em] \displaystyle \int_{\mathbb{R}^d} \pmb{\omega} K(\omega)\mbox{d}\pmb{\omega} & = & \pmb{0}_d, \\[0.8em] \displaystyle \int_{\mathbb{R}^d} \pmb{\omega}\pmb{\omega}^\top K(\omega)\mbox{d}\pmb{\omega} & = & \pmb{I}_d\cdot \end{array} \]

    Se \(K\) for de fato uma densidade de probabilidade multivariada, então as duas últimas equações resumem muitas suposições sobre os kernels marginais, \(\{K_i (\omega_i), \, i = 1,\cdots, d\}\). A segunda equação afirma que as médias dos kernels marginais são todas zero. A terceira equação afirma que os kernels marginais são todos não correlacionados aos pares e que cada um tem variância unitária. Assim, assume-se que qualquer transformação linear simples seja capturada inteiramente na matriz \(H\) e não no kernel.

    Em notação matricial, é possível calcular o erro do estimador de kernel multivariado. Definindo \(\pmb{\omega}=H^{-1}(\pmb{x}-\pmb{y})\), \[ \tag{6.47} \begin{array}{rcl} \mbox{E}\big(\widehat{f}(\pmb{x}) \big) & = & \displaystyle \int_{\mathbb{R}^d} K\big(H^{-1}(\pmb{x}-\pmb{y}) \big) \dfrac{f(\pmb{y})}{|H|}\mbox{d}\pmb{y} \\[0.8em] & = & \displaystyle \int_{\mathbb{R}^d} K\big(\pmb{\omega} \big) f\big(\pmb{x}-H\pmb{\omega}\big)\mbox{d}\pmb{\omega} \\[0.8em] & = & \displaystyle \int_{\mathbb{R}^d} K\big(\pmb{\omega}\big) \left( f(\pmb{x}) -\pmb{\omega}^\top H \nabla f(\pmb{\omega}) +\dfrac{1}{2}\pmb{\omega}^\top H^\top \nabla^2 f(\pmb{\omega})H\pmb{\omega} \right) \mbox{d}\pmb{\omega} \end{array} \] até a segunda ordem.

    Uma simplificação adicional é possível usando a seguinte propriedade do traço (\(\mbox{tr}\)) de uma matriz: \(\mbox{tr}(AB) = \mbox{tr}(BA)\), assumindo que as matrizes \(A\) e \(B\) têm dimensões \(r\times s\) e \(s\times r\), respectivamente. Agora, a forma quadrática na equação (6.47) é uma matriz \(1\times 1\), que trivialmente é igual ao seu traço. Portanto, usando a identidade do traço e trocando as operações de traço e integral, obtemos \[ \mbox{E}\big(\widehat{f}(\pmb{x}) \big) = f(\pmb{x})-0+\dfrac{1}{2}\mbox{tr}\left(\int_{\mathbb{R}^d} \omega\omega^\top K(\pmb{\omega}) \mbox{d}\omega \times H^\top \nabla^2 f(\pmb{x})H \right)\cdot \]

    Como a matriz de covariância de \(K\) é \(\pmb{I}_d\) pela suposição (6.46), o fator integral no traço desaparece. Portanto, \[ \mbox{Viés}\big(\widehat{f}(\pmb{x}) \big) = \dfrac{1}{2}\mbox{tr}\big( H^\top \nabla^2 f(\pmb{x})H\big)=\dfrac{1}{2}\mbox{tr}\big( HH^\top \nabla^2 f(\pmb{x})\big)\cdot \] A seguir, defina o escalar \(h > 0\) e a matriz \(d\times d\) \(A\) satisfazendo \[ H=h\, A \qquad \mbox{onde} \qquad |A|=1\cdot \] Escolher \(A\) para ter determinante igual a 1 significa que a forma elíptica do kernel é inteiramente controlada pela matriz \(AA^\top\) e o tamanho do kernel é inteiramente controlado pelo escalar \(h\).

    Observe que esta parametrização é totalmente geral e permite diferentes parâmetros de suavização para cada dimensão. Por exemplo, se \[ H=\begin{pmatrix} h_1 & \cdots & 0 \\ & \ddots & \\ 0 & \cdots & h_d \end{pmatrix}, \qquad \mbox{temos que} \qquad H=h\times \begin{pmatrix} h_1/h & \cdots & 0 \\ & \ddots & \\ 0 & \cdots & h_d/h \end{pmatrix}, \] onde \(h=(h_1 \, h_2 \cdots h_d)^{1/d}\) é a média geométrica dos \(d\) parâmetros de suavização. Verifique que \(|A|=1\).

    Segue-se que \[ \tag{6.48} \mbox{Viés}\big(\widehat{f}(\pmb{x}) \big) = \dfrac{1}{2} h^2 \mbox{tr}\big( A A^\top \nabla^2 f(\pmb{x})\big), \] para que \[ \mbox{AISB} = \dfrac{1}{4}h^4 \int_{\mathbb{R}^d} \left(\mbox{tr}\big(A A^\top \nabla^2 f(\pmb{x}) \big) \right)^2 \mbox{d}\pmb{x}\cdot \] Como de costume, o termo de variância é dominado por \(\mbox{E}\big(K_H(\pmb{x}-\pmb{x}_i)\big)^2\); portanto, \[ \tag{6.49} \mbox{Var}\big(\widehat{f}(\pmb{x}) \big) = \dfrac{f(\pmb{x})}{n|H|}\int_{\mathbb{R}^d} K(\pmb{\omega})^2 \mbox{d}\pmb{\omega}, \qquad \mbox{o que implica} \qquad \mbox{AIV}=\dfrac{R(K)}{nh^d}\cdot \]


    Teorema 6.5:

    Para um estimador de kernel multivariado geral (6.45) parametrizado por \(H=h A\), \[ \tag{6.50} \mbox{AMISE}=\dfrac{R(K)}{nh^d} +\dfrac{1}{4}h^4 \int_{\mathbb{R}^d} \left(\mbox{tr}\big(A A^\top \nabla^2 f(\pmb{x}) \big) \right)^2 \mbox{d}\pmb{x}\cdot \]


    Demonstração. Scott (2015).


    Apesar das aparências, isso não está usando a mesma largura de banda em cada dimensão, mas sim aplicando um núcleo geral de formato elíptico em uma rotação arbitrária.

    A integral na equação (6.50) será bastante complicada de avaliar, a menos que a matriz \(A\) tenha uma estrutura muito simples. No entanto, Wand and Jones (1995) fornecem uma expressão inteligente na Secção 4.3 que torna muito mais fácil de avaliar esta integral.

    Defina a matriz simétrica \[ M=2\nabla^2 f(\pmb{\pmb{x}})-\mbox{Diag}\big(\nabla^2 f(\pmb{x}) \big)\cdot \] Coloque a porção triangular inferior de \(M\) em um vetor \(\pmb{n}\) de comprimento \(d(d+1)/2\), um procedimento conhecido como operação de meia vetorização, \(\pmb{n} = \mbox{vech}(M)\).

    Sua expressão final é o escalar \[ \tag{6.51} \big( \mbox{vech}(MM^\top)\big)^\top \left( \int \pmb{n}\pmb{n}^\top \mbox{d}\pmb{x}\right) \mbox{vech}(MM^\top)\cdot \] A integral é aplicada aos elementos da matriz \(\pmb{n}\pmb{n}^\top\) de dimensão \(d(d+1)/2 \times d(d+1)/2\).

    As integrais podem ser calculadas ou estimadas usando integração por partes em muitos casos. Como \(AA^\top\) também é uma matriz \(d\times d\) simétrica, \(\mbox{vech}(MM^\top)\) é um vetor de comprimento \(d(d+1)/2\). Veja Wand and Jones (1995) e Exercício 38 para mais detalhes.


    6.3.3 Núcleos de fronteira para regiões irregulares


    Staniswalis, Messer, and Finston (1993) mostraram que um núcleo de fronteira para um domínio arbitrariamente complexo pode ser construído por um dispositivo simples. Suponha que se deseje uma estimativa em \(\pmb{x}\). Eles propõem o uso de um núcleo com suporte esférico, de raio \(h\).

    Apenas as amostras \(\pmb{x}_i\) na esfera de raio \(h\) em torno de \(\pmb{x}\) influenciam a estimativa. Para cada amostra \(\pmb{x}_i\) nessa região, determinar se o diâmetro no qual ela se encontra, sendo o centro o ponto de estimativa \(\pmb{x}\), intersecta a fronteira. Se sim, construa um núcleo de fronteira unidimensional ao longo desse diâmetro. Repetindo essa construção para todas as amostras na esfera, os autores provam que a estimativa resultante mantém a ordem correta do viés.


    6.4 Generalidade do método kernel


    6.4.1 Método delta


    Walter and Blum (1979b) catalogaram a característica comum da já crescente lista de diferentes estimadores de densidade. Ou seja, cada um poderia ser reexpresso como um estimador de kernel. Tal demonstração para estimadores de séries ortogonais foi dada na Seção 6.1.3.

    Reexaminando a equação (3.2), até mesmo o histograma pode ser pensado como um estimador de kernel. Surpreendentemente, este resultado mostrou-se válido mesmo para estimadores que eram soluções para problemas de otimização. Por exemplo, considere um dos vários critérios de máxima verossimilhança penalizada (MPL) sugeridos por Good and Gaskins (1972): \[ \tag{6.52} \max_f \left(\sum_{i=1}^n \log\big(f(x_i\big) -\alpha \int_{-\infty}^\infty f'(x)^2 \mbox{d}x\right) \qquad \mbox{para algum} \qquad \alpha>0\cdot \]

    Sem o termo de penalização de rugosidade em (6.52), a solução seria a função de densidade de probabilidade empírica. Montricher, Tapia, and Thompson (1975) e Klonias (1982) demonstraram que muitos estimadores MPL são estimadores de kernel. A forma das soluções de kernel difere no fato de que os pesos nos kernels não são todos iguais a \(1/n\).

    Para alguns outros algoritmos de estimação de densidade, o kernel equivalente tem peso \(1/n\), mas possui uma largura de banda adaptativa. Um exemplo simples desse tipo é o estimador de \(k\)-ésimo vizinho mais próximo (\(k\)-NN). A estimativa \(k\)-NN em \(x\) é equivalente a uma estimativa de histograma com um intervalo centrado em \(x\), com largura de intervalo suficientemente grande para que o intervalo contenha \(k\) pontos, em duas e três dimensões, a forma do intervalo do histograma é um círculo e uma esfera, respectivamente. Assim, o kernel equivalente em \(\mathbb{R}^d\) é simplesmente uma densidade uniforme sobre a bola unitária em \(\mathbb{R}^d\), mas com larguras de intervalo que se adaptam a \(x\).


    6.4.2 Teorema do kernel generalizado


    Existe uma base teórica para as observações empíricas de Walter e Blum de que muitos algoritmos podem ser vistos como estimativas de kernel generalizadas. Terrell apresentou um teorema nesse sentido que contém um algoritmo construtivo para obter o kernel generalizado de qualquer estimador de densidade (ver G. R. Terrell and Scott 1992a).

    A construção não se limita a estimadores não paramétricos, um fato que será explorado posteriormente.


    Teorema 6.6:

    Qualquer estimador de densidade que seja um funcional contínuo e diferenciável segundo Gâteaux na função de distribuição empírica pode ser escrito como \[ \tag{6.53} \widehat{f}(x)=\dfrac{1}{n}\sum_{i=1}^n K(x,x_i,F_n), \] onde \(K\) é a derivada Gâteaux de \(\widehat{f}\) sob variação de \(x_i\).


    Demonstração. Scott (2015).


    A derivada de Gâteaux de um funcional \(T\) na função \(\phi\) na direção da função \(\eta\) é definida como sendo \[ \tag{6.54} \mbox{DT}(\phi)[\eta]=\lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \Big(T(\phi+\epsilon\eta)-T(\phi) \Big)\cdot \]

    O Teorema 6.6, que é demonstrado abaixo, possui uma versão multivariada análoga (G. R. Terrell and Scott 1992a). O núcleo \(K\) simplesmente mede a influência de \(x_i\) em \(\widehat{f}(x)\). À medida que \(F_n\) converge para \(F\), então, assintoticamente, a forma de \(K\) é independente das \(n-1\) observações restantes.

    Assim, qualquer estimador de densidade contínua pode ser escrito (assintoticamente) como \[ \tag{6.55} \widehat{f}(x)=\dfrac{1}{n}\sum_{i=1}^n K(x,x_i,\eta)=\dfrac{1}{n}\sum_{i=1}^n K_n(x,x_i)\cdot \]


    6.4.2.1 Demonstração do resultado geral do kernel


    A função de distribuição acumulada empírica na equação (2.1) pode ser escrita na forma incomum \[ \tag{6.56} F_n(\cdot)=\dfrac{1}{n}\sum_{i=1}^n \pmb{I}_{[x_i,\infty)}(\cdot)\cdot \]

    Escrevemos o estimador de densidade como um operador \(\widehat{f}(x)=T_x\big(F_n\big)\). Definindo \[ \tag{6.57} \begin{array}{rcl} K(x,y,F_n) & = & \displaystyle \lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \left( T_x\big( (1-\epsilon)F_n+\epsilon\pmb{I}_{[y,\infty)}\big) - (1-\epsilon)T_x\big(F_n \big) \right) \\[0.8em] & = & \displaystyle \lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \left( T_x\Big( F_n+\epsilon\big(\pmb{I}_{[y,\infty)}-F_n\big)\Big) - T_x\big(F_n \big) \right)+T_x\big(F_n \big) \\[0.8em] & = & \mbox{DT}(F_n)\left(\pmb{I}_{[y,\infty)}-F_n \right)+\widehat{f}(x), \end{array} \] onde \(\mbox{DT}(\phi)[\eta]\) é a derivada de Gâteaux de \(T\) em \(\phi\) na direção \(\eta\). A proposição (2.7) de Tapia (1971) mostra que a derivada de Gâteaux é linear em seu segundo argumento, então \[ \begin{array}{rcl} \displaystyle \dfrac{1}{n} \sum_{i=1}^n K(x,x_i,F_n) & = & \displaystyle \dfrac{1}{n}\sum_{i=1}^n \mbox{DT}_x (F_n)\left( \pmb{I}_{[x_i,\infty)}-F_n \right) +\widehat{f}(x) \\[0.8em] & = & \displaystyle \mbox{DT}_x(F_n)\left(\dfrac{1}{n}\sum_{i=1}^n \pmb{I}_{[x_i,\infty)} -F_n \right)+\widehat{f}(x)\\[0.8em] & = & 0+\widehat{f}(x), \end{array} \] onde o termo entre parênteses é 0 pela equação (6.56).

    Observe que, por linearidade, a variação de Gâteaux na direção 0, \(\mbox{DT}(\phi)[0]\), é identicamente 0. Isso conclui a demonstração.


    6.4.2.2 Caracterização de um estimador não paramétrico


    Um estimador é definido como não paramétrico quando é consistente na média quadrática para uma ampla classe de funções de densidade.

    Com um pouco de esforço, essa definição se traduz em requisitos específicos para o kernel equivalente, \(K_n(x,y)\). A partir dos muitos exemplos anteriores, um estimador não paramétrico que seja consistente deve ser local; isto é, a influência de pontos amostrais fora de uma vizinhança de \(x\) deve desaparecer quando \(n\to\infty\).

    Como sugerido pelo termo “sequências de funções delta”, cunhado por Watson and Leadbetter (1963), \[ \lim_{n\to\infty} K_n(x,x_i,F_n)=\delta(x-x_i)\cdot \] No entanto, \(K_n\) não deve ser igual a \(\delta(x-x_i)\) para \(n\) finito, como foi o caso na Seção 6.1.3 com o estimador de série ortogonal com todos os coeficientes de \(\widehat{f}_\nu\) incluídos.

    Começando com o viés da forma do kernel assintótico em (6.55) \[ \begin{array}{rcl} \mbox{E}\big(\widehat{f}\big) & = & \displaystyle \mbox{E}\big(K_n(x,X) \big) = \int K_n(x,y)f(y)\mbox{d}y \\[0.8em] & = & \displaystyle \int K_n(x,y) \left( f(x)+(y-x)f'(x)+\dfrac{(y-x)^2}{2} f''(\xi_y)\right) \mbox{d}y, \end{array} \] pela versão exata do teorema de Taylor, onde \(\xi_y\in (x,y)\).

    Para ser assintoticamente não viesado, todos os três termos a seguir devem ser nulos: \[ \tag{6.58} \begin{array}{rcl} \mbox{Viés}\big(\widehat{f}\big) & = & \displaystyle f(x)\left( \int K_n(x,y)\mbox{d}y-1\right) + f'(x)\int K_n(x,y)(y-x)\mbox{d}y \\[0.8em] & & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \displaystyle \dfrac{1}{2}\int K_n(x,y)(y-x)^2f''(\xi_y)\mbox{d}y\cdot \end{array} \] Portanto, a primeira condição é que \[ \lim_{n\to\infty} \int K_n(x,y)\mbox{d}y=1 \qquad \forall x\in \mathbb{R}^1\cdot \] Suponha que o estimador tenha sido reescalado de forma que a integral seja exatamente 1 para todo \(n\). Portanto, defina a variável aleatória \(Y\) com função densidade de probabilidade \(K_n(x,\cdot)\). A segunda condição é que \[ \lim_{n\to\infty} \int K_n(x,y)y\mbox{d}y = \int K_n(x,y)x\mbox{d}x = x \qquad \forall x, \] o que implica \[ \lim_{n\to\infty} K_n(x,y)=\delta(y-x), \] como sugeriram Watson and Leadbetter (1963) e Walter and Blum (1979a).

    O comportamento preciso do viés é determinado pela taxa na qual isso ocorre. Por exemplo, suponha que \[ \int K_n(x, y) y\mbox{d}y = x \] para todo \(n\) e que \[ \tag{6.59} \sigma_{x,n}^2 =\int K_n(x,y)(y-x)^2 \mbox{d}y \neq 0 \] para \(n\) finito, de modo que os dois primeiros momentos das variáveis aleatórias \(Y\sim (x,\sigma_{x,n}^2)\) e \(T=(Y-x)/\sigma_{x,n}\sim (0,1)\).

    Suponha que a função densidade de \(T\), que é uma transformação linear simples de \(K_n\), convirja para uma densidade adequada: \[ \widetilde{L}_n(x,t)=K_n(x,x+\sigma_{x,n}t)\sigma_{x,n} \longrightarrow L(x,t) \qquad \mbox{quando} \qquad n\to\infty\cdot \] Então, o último termo de viés em (6.58) pode ser aproximado por \[ \begin{array}{rcl} \displaystyle \dfrac{1}{2}f''(x) \int K_n(x,y)(y-x)^2 \mbox{d}y & = & \displaystyle \dfrac{1}{2}f''(x)\int K_n(x,x+\sigma_{x,n} t) (\sigma_{x,n} t)^2 \mbox{d}t \\[0.8em] & = & \displaystyle \dfrac{1}{2}\sigma_{x,n}^2 f''(x)\int t^2 \widetilde{L}(x,t)\mbox{d}t \\[0.8em] & \approx & \displaystyle \dfrac{1}{2}\sigma_{x,n}^2 f''(x)\int t^2 L(x,t)\mbox{d}t = \dfrac{1}{2}\sigma_{x,n}^2 f''(x) \\[0.8em] \end{array} \] como \(\mbox{Var}(T) = 1\), o viés é \(O(\sigma_{x,n}^2)\), a taxa usual para kernels de segunda ordem. Portanto, a terceira condição é que \(\sigma_{x,n}\to 0\) quando \(n\to\infty\).

    Para que a variância se anule assintoticamente, considere \[ \begin{array}{rcl} \mbox{Var}\big(\widehat{f}\big) & = & \dfrac{1}{n}\mbox{Var}\big(\widehat{f}\big) \leq \mbox{E}\big(K_n(x,X)^2\big) = \dfrac{1}{n} \displaystyle \int K_n(x,y) f(y)\mbox{d}y \\[0.8em] & = & \displaystyle \dfrac{1}{n}\int K_n(x,x+\sigma_{x,n}t)^2 f(x+\sigma_{x,n}t)\sigma_{x,n} \mbox{d}t \\[0.8em] & = & \displaystyle \dfrac{1}{n\, \sigma_{x,n}} \int \widetilde{L}_n(x,t)^2 \big(f(x)+\cdots \big) \mbox{d}t \approx \dfrac{f(x) R\big(L(x,\cdot) \big)}{n\, \sigma_{x,n}}\cdot \end{array} \]

    Assim, a quarta condição exigida é que a variância do kernel equivalente satisfaça \(n\, \sigma_{x,n}\to\infty\) quando \(n\to\infty\).


    Teorema 6.7:

    Suponha que \(\widehat{f}\) seja um estimador de densidade com núcleo assintótico equivalente \(K_n (x,y)\) e que \(\sigma_{x,n}^2\) definido em (6.59) seja limitado e não nulo. Então \(\widehat{f}\) é um estimador de densidade não paramétrico se, para todo \(x\in\mathbb{R}^1\), \[ \tag{6.60} \begin{array}{rcl} \displaystyle \lim_{n\to\infty} \int K_n(x,y)\mbox{d}y & = & 1, \\[0.8em] \displaystyle \lim_{n\to\infty} \int K_n(x,y)\, y\, \mbox{d}y & = & x, \\[0.8em] \displaystyle \lim_{n\to\infty} \sigma_{x,n} & = & 0, \\[0.8em] \displaystyle \lim_{n\to\infty} n\, \sigma_{x,n} & = & 1\cdot \end{array} \]


    Demonstração G. R. Terrell (1984).

    6.4.2.3 Núcleos equivalentes de estimadores paramétricos


    O Teorema 6.6 mostra como construir o núcleo para qualquer estimador de densidade, paramétrico ou não paramétrico. Por exemplo, considere a estimativa paramétrica de \[ f = N(\mu,1) = \phi(x \, |\, \mu,1) \qquad \mbox{com} \qquad \widehat{\mu}=\overline{x} \cdot \] Assim, \(T_x(F_n) = \phi( x \, | \, \overline{x}, 1)\). Examinando o argumento na primeira linha da equação (6.57) e comparando-o com a definição da função de distribuição empírica (ecdf) em (6.56), torna-se evidente que a ecdf \[ \dfrac{1}{n} \sum_{i=1}^n \delta (x-x_i) \] está sendo substituída por \[ \dfrac{1-\epsilon}{n}\sum_{i=1}^n \delta(x-x_i)+\epsilon \, \delta(x-y) \] que é a função de probabilidade empírica original com uma pequena porção da massa de probabilidade removida proporcionalmente e colocada em \(x = y\). A média amostral dessa função de densidade de probabilidade empírica perturbada é \((1-\epsilon)\overline{x}+\epsilon y\).

    Assim, o kernel pode ser calculado diretamente a partir da equação (6.57) por \[ \tag{6.61} \begin{array}{rcl} K(x,y,F_n) & = & \displaystyle \lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \left( \phi\big(x \, | \, (1-\epsilon)\overline{x}+\epsilon \, y,1 \big) -(1-\epsilon)\phi(x \, | \, \overline{x},1)\right) \\[0.8em] & = & \displaystyle \dfrac{1+(y-\overline{x})(x-\overline{x})}{\sqrt{2\pi}} \exp\left( -\frac{1}{2}(x-\overline{x})^2\right) \end{array} \] por uma série de Taylor (ver Erercício 28).

    O kernel assintótico equivalente é \[ K(x,y)=\lim_{n\to\infty} K(x,y,F_n) + \dfrac{1+(y-x)(x-\mu)}{\sqrt{2\pi}} \exp\left(-\frac{1}{2}(x-\mu)^2 \right)\cdot \]

    Este kernel nunca é local e, portanto, o estimador não é não paramétrico, caso houvesse alguma dúvida. Observe que o estimador de kernel paramétrico com kernel (6.61) é bastante bom, pois \[ \widehat{f}(x)=\dfrac{1}{n}\sum_{i=1}^n \dfrac{1+(x_i-\overline{x})(x-\overline{x})}{\sqrt{2\pi}}\exp\left( -\frac{1}{2}(x-\overline{x})^2\right)=\dfrac{1}{\sqrt{2\pi}}\exp\left( -\frac{1}{2}(x-\overline{x})^2\right)\cdot \]

    Naturalmente, tal como no contexto paramétrico, este estimador de “kernel” será sempre \(\phi(x \, | \, \overline{x}, 1)\) para todos os conjuntos de dados, independentemente de quão não normal seja a densidade subjacente.


    6.5 Validação cruzada


    6.5.1 Dados univariados


    O objetivo é ir além dos resultados teóricos para a especificação de largura de banda ideal e alcançar algoritmos práticos baseados em dados. As regras de referência supersuavizadas e normais fornecem escolhas iniciais razoáveis, juntamente com as modificações simples baseadas na assimetria e curtose da amostra, apresentadas na Seção 4.1.2.

    Os fatores de modificação da largura de banda para kernels positivos e o polígono de frequência são idênticos. Os algoritmos de validação cruzada não viesada e viesada para o histograma são facilmente estendidos tanto para estimadores de kernel quanto para estimadores ASH. No entanto, o desenvolvimento de algoritmos de seleção de largura de banda para estimadores de kernel avançou além daqueles para o histograma.

    A importância da escolha da largura de banda ideal é facilmente superestimada. A partir dos resultados de sensibilidade na Tabela 3.3, qualquer escolha de \(h\) dentro de 15-20% de \(h^*\) geralmente será suficiente. G. R. Terrell (1990) chegou a sugerir que uma regra de suavização excessiva fosse geralmente usada.

    Com dados reais, é relativamente fácil examinar uma sequência de estimativas com base na sequência de parâmetros de suavização \[ h=\widehat{h}_{OS}/1.05^k \qquad k=0,1,2,\cdots, \] começando com a largura de banda de suavização excessiva amostral \(\widehat{h}_{OS}\) e parando quando a estimativa apresentar alguma instabilidade e ruído muito local próximo aos picos.

    Silverman (1978a) caracterizou a quantidade esperada de ruído local presente em \(\widehat{f}''(x)\) na norma \(L_\infty\), ver equação (6.65). Ele sugeriu examinar gráficos de \(\widehat{f}''(x)\) interativamente, além de \(\widehat{f}(x)\), um procedimento conhecido como método do gráfico de teste.

    O algoritmo de validação cruzada enviesada, que estima \(R(f'')\) a partir de uma estimativa de kernel, tenta corrigir a superestimação por \(R(\widehat{f}''(\cdot))\) que é o resultado direto da presença do ruído local, ver Seção 6.5.1.3. No entanto, o poder da abordagem interativa para a seleção de largura de banda não deve ser subestimado. O processo interativo pode ser bastante simples em sistemas que suportam animação, como LISP-STAT (Tierney 1990).

    Deve-se enfatizar que, de um ponto de vista exploratório, todas as escolhas da largura de banda \(h\) levam a estimativas de densidade úteis. Larguras de banda grandes, por um lado, fornecem uma imagem da estrutura global na densidade desconhecida, incluindo características gerais como assimetria, outliers, clusters e localização. Larguras de banda pequenas, por outro lado, revelam estruturas locais que podem ou não estar presentes na densidade real.

    Além disso, a otimalidade de \(h\) depende não apenas da escolha da métrica \(L_p\), mas também da característica na densidade a ser enfatizada \(F\), \(f\), \(f'\), \(f''\), \(\cdots\). No entanto, o desejo por procedimentos de seleção de largura de banda totalmente automáticos e confiáveis levou inevitavelmente a uma série de novos algoritmos. Esses procedimentos objetivos (não subjetivos) geralmente têm o objetivo declarado de encontrar uma largura de banda que minimize o erro \(L_2\) real, em vez de usar a largura de banda que minimize o erro \(L_2\) esperado (\(\mbox{MISE}\)).

    Em um artigo inicial, Wahba (1981) expressou a expectativa de que seu algoritmo de validação cruzada generalizado atingiria esse objetivo. A melhor largura de banda \(L_2\), \(h_{ISE}\), permaneceu o alvo na validação cruzada não enviesada para Hall and Marron (1987b), Hall and Marron (1987a). Scott and Factor (1981) expressaram a opinião de que \(h_{MISE}\) era um alvo apropriado. A largura de banda ótima para \(\mbox{MISE}\) depende apenas de \((f_n)\), enquanto a largura de banda ótima para \(\mbox{ISE}\) também depende da amostra, \((f, x, \{x_i\})\). No entanto, foi demonstrado que a correlação da amostra entre \(h_{ISE}\) e se aproximou de -0.70 para dados normais (Scott and Terrell 1987; Scott 1988). Dada uma correlação negativa tão grande com a escala dos dados, rastrear \(h_{ISE}\) de perto exigiria adivinhar se \(\widehat{\sigma}>\sigma\), ou vice-versa, uma tarefa difícil.

    Um fato interessante é que, embora \(\widehat{h}_{ISE}\approx h_{MISE}\), \[ \dfrac{\widehat{\sigma}_{\widehat{h}_{ISE}}}{h_{MISE}} = O(n^{-1/10}), \] de modo que a largura de banda ótima para o \(\mbox{ISE}\) converge lentamente para o \(h_{MISE}\), como demonstrado por Hall and Marron (1987b). Scott and Terrell (1987) mostraram que as larguras de banda da validação cruzada não enviesada (\(\mbox{UCV}\)) e da validação cruzada enviesada (\(\mbox{BCV}\)) convergiam para o \(h_{MISE}\) na mesma taxa lenta.

    Algumas das extensões mais recentes conseguiram elevar a taxa de convergência relativa até \(O(n^{-1/2})\), que é a melhor taxa possível (Hall and Marron 1991). Esses procedimentos requerem a introdução de um ou dois parâmetros de suavização auxiliares. Dada a lenta taxa de convergência do inatingível \(h_{ISE}\), talvez não esteja claro se existe alguma vantagem prática em taxas mais rápidas. Essa questão é examinada mais detalhadamente na Seção 6.5.1.5.

    Em um estudo empírico do desempenho de nove algoritmos de validação cruzada e quatro densidades de amostragem, Jones and Kappenman (1992) relataram a “ampla equivalência de quase todos” esses algoritmos em relação ao \(\mbox{ISE}\) observado. Outras simulações (Scott and Factor 1981; Bowman 1985; Scott and Terrell 1987; Park and Marron 1990) relataram menor equivalência entre os próprios valores estimados dos parâmetros de suavização.

    Jones and Kappenman (1992) relataram que a largura de banda \(\mbox{AMISE}\) fixa \(h^*\) superou todos os nove algoritmos de validação cruzada em relação ao \(\mbox{ISE}\). Seus resultados reforçam a sugestão de que \(h^*\) é um alvo apropriado e que qualquer escolha dentro de 15-20% de \(h^*\) deve ser adequada. A maioria dos esforços agora se concentra no \(h_{MISE}\) como a largura de banda alvo.


    6.5.1.1 Primeiros esforços na seleção de largura de banda


    As primeiras ideias de seleção de largura de banda baseadas em dados ou amostrais surgiram no contexto de estimadores de séries ortogonais, que foram discutidos na Seção 6.1.3.

    Usando os pesos de Tarter-Kronmal (6.8) e a representação na equação (6.5), o erro pontual é \[ \widehat{f}(x)-f(x)=\sum_{\nu=-m}^m \widehat{f}_\nu \phi_\nu(x) -\sum_{\nu=-\infty}^\infty f_\nu\phi_\nu(x)\cdot \]

    Como as funções de base são ortonormais, \[ \mbox{ISE}=\sum_{\nu=-m}^m || \widehat{f}_\nu -f_\nu ||^2 + \sum_{\nu\notin [-m,m]} || f_\nu ||^2 \cdot \]

    Lembre-se que \(\mbox{MISE} = \mbox{E}\big(\mbox{ISE}\big)\). O procedimento de seleção de Tarter e Kronmal forneceu estimativas não viesadas do incremento no \(\mbox{MISE}\) ao passar de uma série com \(m-1\) termos para uma com \(m\) termos, observando a igualdade dos termos \(\pm \nu\) do \(\mbox{MISE}\): \[ \tag{6.62} \mbox{MISE}(m)-\mbox{MISE}(m-1) = 2 \left( \mbox{E}\big( ||\widehat{f}_m -f_m ||^2 \big) - || f_m||^2\right)\cdot \]

    Estimativas não viesadas dos dois termos do lado direito podem ser obtidas para o estimador de Fourier na equação (6.6), como \(\mbox{E}\big(\widehat{f}_\nu\big) = f_\nu\) e \[ \mbox{E}\big(\widehat{f}_\nu\widehat{f}_\nu^* \big) = \dfrac{1-(n-1)| \widehat{f}_\nu|^2}{n}, \] onde \(\widehat{f}_\nu^*\) denota o conjugado complexo de \(\widehat{f}\nu\); portanto, \[ \mbox{Var}\big(\widehat{f}_\nu\big) = \mbox{E}\big(\widehat{f}_\nu\widehat{f}_\nu^* \big)-| \widehat{f}_\nu|^2 = \dfrac{1-| \widehat{f}_\nu|^2}{n} \cdot \]

    A escolha baseada em dados para \(m\) é alcançada quando o incremento se torna positivo. Observe que aceitar a inclusão do \(m\)-ésimo coeficiente no estimador da série resulta do julgamento de que a variância adicional de \(\widehat{f}_m\) é menor que a redução no viés \(|| f_m||^2\).

    Normalmente, menos de seis termos são escolhidos, de modo que apenas ajustes relativamente grosseiros possam ser feitos na suavidade do estimador de densidade. Às vezes, o algoritmo omite termos de ordem superior. Mas a verdadeira importância deste algoritmo reside na sua afirmação de ser o primeiro algoritmo de validação cruzada não enviesado para um estimador de densidade.

    Da mesma forma, o crédito pelo primeiro algoritmo de validação cruzada enviesado é atribuído a Wahba (1981) com seu algoritmo de validação cruzada generalizado. Ela usou as mesmas estimativas não enviesadas dos coeficientes de Fourier que Tarter e Kronmal, mas com sua escolha mais suave de pesos, perdeu a simplicidade de examinar as mudanças incrementais no \(\mbox{MISE}\). No entanto, essas mesmas estimativas não enviesadas dos coeficientes de Fourier levam a uma boa estimativa do \(\mbox{AMISE}\). Ao ignorar todos os coeficientes de Fourier desconhecidos para \(|\nu| > n/2\), um pequeno viés é introduzido. Ambos os grupos recomendam mostrar graficamente a função de risco estimada para encontrar o melhor parâmetro de suavização baseado em dados, em vez de recorrer à otimização numérica (cega).

    A primeira tentativa de escolher o parâmetro de suavização do kernel de forma totalmente automática foi um algoritmo de máxima verossimilhança modificado, proposto por Habbema, Hermans, and Van Der Broek (1974) e Duin (1976). Embora não tenha resistido ao teste do tempo, é significativo por ter introduzido uma modificação do tipo “deixar um de fora” ao critério usual de máxima verossimilhança (ML).

    A escolha da largura de banda \(h\) para maximizar o critério usual de ML resulta na seguinte função de densidade de probabilidade (pdf) empírica (aproximada): \[ \arg \max_h \sum_{i=1}^n \log\big(\widehat{f}(x_i,h)\big) =0 \qquad \mbox{que implica} \qquad \widehat{f}(x_i,h=0)=\dfrac{1}{n}\sum_{i=1}^n \delta(x-x_i), \] uma solução com versossimilhança “infinita”.

    O problema surge porque, quando \(h\to 0\), a contribuição do próprio ponto \(x_i\) para a verossimilhança em \(x = x_i\) torna-se infinita. Os autores procuraram eliminar essa “autocontribuição” modificando o critério de máxima verossimilhança: \[ \max_h \sum_{i=1}^n \log\big( \widehat{f}_{-i}(x_i;h)\big), \] onde \(\widehat{f}_{-i}(x_i;h)\) é um estimador de kernel baseado nos \(n-1\) pontos de dados, excluindo \(x_i\), e então avaliado nesses pontos.

    Apesar de alguns resultados empíricos promissores em amostras pequenas e de consistência (Chow, Geman, and Wu 1983), o algoritmo mostrou-se excessivamente influenciado por outliers e clusters compactos (Schuster and Gregory 1981; Scott and Factor 1981).

    Com um kernel de suporte finito, por exemplo, a largura de banda não pode ser menor que \(x_{(2)}-x_{(1)}\), que é a distância entre as duas primeiras estatísticas de ordem; para muitas densidades, a distância entre essas estatísticas de ordem não converge para zero e, portanto, a largura de banda não converge para zero como requerido.

    Um algoritmo de ponto fixo simples foi sugerido por Scott, Tapia, and Thompson (1977). Para estimativas de kernel, a única incógnita em \(h^*\) em (6.18) é \(R(f'')\). Se um kernel normal for usado, um cálculo direto mostra que \[ \tag{6.63} R\big(\widehat{f}''_h\big) = \dfrac{3}{8\sqrt{\pi} \, n^2 h^5} \sum_{i=1}^n \sum_{j=1}^n \left( -\Delta_{ij}^2 +\dfrac{1}{12} \Delta_{ij}^4\right) \exp\Big(-\frac{1}{4}\Delta_{ij}^2\Big), \] onde \(\Delta_{ij}=(x_i-x_j)/h\).

    Seguindo a equação (6.18), a busca por um valor de ponto fixo para \(h^*\) é realizada por meio de iterações \[ h_{k+1}=\left( \dfrac{R(K)}{n \, \sigma_K^4 R\big( \widehat{f}_{h_k}''\big)}\right)^{1/5}, \] com \(h_0\) escolhido como a largura de banda de referência normal.

    Como a razão entre as larguras de banda ótimas para estimar \(f\) e \(f''\) diverge quando \(n\to\infty\), fica claro que o algoritmo não é consistente. O fato de o algoritmo ter funcionado bem para amostras pequenas (Scott and Factor 1981) não é surpreendente, visto que as larguras de banda ótimas são razoavelmente próximas umas das outras para amostras pequenas. Observe que este algoritmo, conforme apresentado, não fornece nenhuma estimativa do \(\mbox{MISE}\).

    É simples usar a estimativa de rugosidade (6.63) na equação (6.18), seguindo a sugestão de Wahba, para obter \[ \tag{6.64} \widehat{\mbox{AMISE}}(h) = \dfrac{R(K)}{n\, h}+\dfrac{1}{4}\sigma_K^4 h^4 R\big(\widehat{f}_h'' \big)\cdot \]

    Encontrar o minimizador da equação (6.64) não só fornece uma estimativa da largura de banda amostral, como também uma estimativa do \(\mbox{MISE}\). Essa ideia foi ressuscitada com validação cruzada enviesada usando um estimador consistente para \(R(f'')\) em vez da equação (6.63).

    Alternativamente, uma largura de banda muito maior, apropriada para \(f''\) em vez de \(f\), pode ser usada na iteração. Sheather (1983) propôs um esquema desse tipo ao tentar estimar a densidade na origem, fornecendo um exemplo dos estimadores plug-in (PI) discutidos por Woodroofe (1970). A motivação de Sheather é discutida em Sheather and Jones (1991). Em uma palestra em 1991, Gasser também relatou sucesso dessa maneira na escolha de uma largura de banda global, inflando a largura de banda usada em \(R\big(\widehat{f}''_h\big)\) pelo fator \(n^{-1/10}\).

    Para um núcleo normal, Silverman (1978a) provou que \[ \tag{6.65} \dfrac{\sup \left| \widehat{f}''- \mbox{E}\big(\widehat{f}'' \big) \right| }{\sup \left| \mbox{E}\big(\widehat{f}'' \big)\right| } \approx 0.4\cdot \]

    Ele propôs escolher a largura de banda onde a relação entre ruído e sinal parece ser de 0.4. Essa relação é diferente para outros kernels. O procedimento do gráfico de teste também pode ser usado no caso bivariado.


    6.5.1.2 Suavização excessiva


    A derivação da regra de suavização excessiva para estimadores de kernel será construtiva, diferentemente das derivações anteriores das regras de suavização excessiva para histograma e polígono de frequência. A versão preferida utiliza a variância como medida de escala. Outras medidas de escala foram consideradas por Terrell G. R. and Scott (1985) e G. R. Terrell (1990).

    Considere o problema variacional \[ \tag{6.66} \min_f \int_{-\infty}^\infty f''(x)^2 \mbox{d}x \qquad \mbox{restrito à} \qquad \int_{-\infty}^\infty f(x)\mbox{d}x=1 \; \mbox{e} \; \int_{-\infty}^\infty x² f(x)\mbox{d}x=1\cdot \] Claramente, a solução será simétrica. O Lagrangiano associada é \[ L(f)=\int_{-\infty}^\infty f''(x)^2 \mbox{d}x+\lambda_1 \left(\int_{-\infty}^\infty f(x)\mbox{d}x-1 \right)+\lambda_2 \left( \int_{-\infty}^\infty x^2 f(x)\mbox{d}x-1\right)\cdot \]

    Como uma solução, a variação de Gâteaux, definida na equação (6.54), do Lagrangiano em qualquer “direção” \(\eta\) deve se anular. Por exemplo, a variação de Gâteaux de \[ \Psi(f)=\int_{-\infty}^\infty f''(x)^2 \mbox{d}x \] é \[ \begin{array}{rcl} \Psi'(f)[\eta] & = & \displaystyle \lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \left( \int_{-\infty}^\infty \Big( \big(f(x)+\epsilon \eta(x) \big)''\Big)^2 \mbox{d}x \right) - \int_{-\infty}^\infty \big( f''(x)\big)^2 \mbox{d}x\\[0.8em] & = & \displaystyle \lim_{\epsilon\to 0} \dfrac{1}{\epsilon} \left( \int_{-\infty}^\infty \Big( \big(f''(x)\big)^2+2\epsilon f''(x) \eta''(x) +\epsilon^2 \big(\eta'' \big)^2 - \big(f''(x) \big)^2\Big) \mbox{d}x \right) = \int_{-\infty}^\infty 2 \big( f''(x)\big)^2\eta''(x) \mbox{d}x\cdot \end{array} \]

    Calculando a variação de Gâteaux de \(L(f)\), temos \[ \tag{6.67} \begin{array}{rcl} 0 & = & \displaystyle L'(f)[\eta] = \int_{-\infty}^\infty 2f''(x) \eta''(x)\mbox{d}x+\lambda_1 \int_{-\infty}^\infty \eta(x)\mbox{d}x +\lambda_2\int_{-\infty}^\infty x^2 \eta(x)\mbox{d}x \\[0.8em] & = & \displaystyle 2f''(x)\eta'(x)\Bigg|_{-c}^c -2f'''(x)\eta(x)\Bigg|_{-c}^c+\int_{-\infty}^\infty \left(2f^{iv}(x)+\lambda_1 +\lambda_2 x^2 \right)\eta(x)\mbox{d}x \end{array} \] após integrar por partes duas vezes e onde \(c\) é o limite, possivelmente infinito, para \(f\).

    Agora, η(±c) deve se anular para que não haja descontinuidade na solução; portanto, o segundo termo se anula. Os dois termos restantes devem se anular para todas as escolhas possíveis de η; portanto, f(±c) deve se anular, restando apenas a integral. Segue-se que o integrando deve se anular e que f é um polinômio de sexta ordem com apenas potências pares (por simetria).

    Portanto, a solução deve assumir a forma \[ f(x)=a(x-c)^3 (x+c)^3 \] de modo que \(f''(\pm c)=0\).

    As duas restrições em (6.66) impõem duas condições lineares às incógnitas \((a,c)\), resultando em que \[ f^*(x)=\dfrac{35}{69.984}(9-x^2)^3_{+} \qquad \mbox{e} \qquad R\big( (f^*)'' \big)=\dfrac{35}{243}\cdot \]

    Uma simples mudança de variáveis mostra que \(R(f''')\geq 35/(243 \, \sigma^5)\); portanto, \[ \tag{6.68} h^* = \left( \dfrac{R(K)}{n\, \sigma_K^4 R(f'')}\right)^{1/5}\leq \left( \dfrac{243 \, \sigma^5 R(K)}{35 \, n \, \sigma_K^4}\right)^{1/5} \] o qual implica em \[ \tag{6.69} \mbox{Regra de suavização excessiva:} \qquad h_{OS}=\dfrac{3}{n^{1/5}}\left(\dfrac{R(K)}{35\, \sigma_K^4} \right) \sigma\cdot \]

    Para o kernel normal \(h_{OS}=1.144 \, \sigma \, n^{-1/5}\). Para o kernle biweight, é exatamente \(h_{OS}=3 \, \sigma \, n^{-1/5}\). A régua é 1.08 vezes mais larga que a regra de referência normal.


    6.5.1.3 Validação cruzada não viesada e viesada


    A apresentação nesta seção dependerá fortemente das referências para certos detalhes. A essência é a mesma da aplicação ao histograma e ao polígono de frequência.

    O fato notável é que a justificativa da validação vruzada não viesada (UCV) para o histograma é inteiramente geral. Assim, a definição na equação (3.52) se aplica no caso de um estimador de kernel. Para o caso de um kernel normal, Rudemo (1982) e Bowman (1984) mostraram que, substituindo \(n\pm 1\) por \(n\) para simplificar, \[ \tag{6.70} \mbox{UCV}(h)=\dfrac{1}{2nh\sqrt{\pi}}+\dfrac{1}{n^2 h\sqrt{\pi}}\sum_{i<j} \left( e^{-\Delta^2_{ij}/4}-\sqrt{8} \, e^{-\Delta^2_{ij}/2}\right), \] que é um caso especial da fórmula geral de Scott and Terrell (1987) \[ \mbox{UCV}(h)=\dfrac{R(K)}{nh}+\dfrac{2}{n^2h}\sum_{i<j} \gamma(\Delta_{ij}), \] onde \[ \gamma(\Delta)=\int K(\omega)K(\omega+\Delta)\mbox{d}\omega\cdot \]

    O algoritmo \(\mbox{BCV}\) decorre do resultado de (Scott and Terrell 1987) de que \[ \mbox{ER}\big( \widehat{f}''_h\big) =R(f'')+\dfrac{R(K)}{n\, h^5}+O(h^2), \] onde \(R(K'')/(n\, h^5)\) é assintoticamente uma constante, representando o ruído fixo, mas finito, que existe na estimação kernel.

    Portanto, \(R\big(\widehat{f}''_h\big)-R(K'')/(n\, h^5)\) é um estimador assintoticamente não viesado para a rugosidade desconhecida \(R(f'')\). Substituindo na equação (6.18), a estimativa do \(\mbox{MISE}\) torna-se \[ \tag{6.71} \mbox{BCV}(h)=\dfrac{R(K)}{nh}+\dfrac{\sigma_K^4}{2n^2h}\sum_{i<j} \phi(\Delta_{ij}), \] onde \[ \phi(\Delta)=\int K''(\omega)K''(\omega+\Delta)\mbox{d}\omega\cdot \]

    A semelhança entre as fórmulas gerais de \(\mbox{UCV}\) e \(\mbox{BCV}\) é notável, considerando suas origens bastante diferentes. No caso de um kernel normal, \[ \tag{6.72} \mbox{BCV}(h)=\dfrac{1}{2nh\sqrt{\pi}}+\dfrac{1}{64n^2h\sqrt{\pi}}\sum_{i<j} \big(\Delta_{ij}^4-12\Delta_{ij}^2+12 \big)e^{-\Delta_{ij}^2/4}\cdot \]

    Como antes, \(\displaystyle \lim_{h\to\infty} \mbox{BCV}(h) = 0\); portanto, \(\widehat{h}_{BCV}\) é considerado o maior minimizador local de \(\mbox{BCV}(h)\) menor ou igual à largura de banda suavizada.

    A teoria assintótica dos critérios \(\mbox{CV}\) é direta, uma vez que se reconhece que a parte estocástica está toda contida nas chamadas \(U\)-estatísticas, que são somas duplas da forma \[ U_n = \sum_{i<j} H_n(X_i,Y_j)\cdot \]

    Hall (1984) provou que se a função \(H_n\) for simétrica e a variável aleatória \(\mbox{E}\big(H_n(X,Y)\, |\, X\big) = 0\), então, juntamente com uma certa condição de momento, \[ U_n=\mbox{AN}\left(0,\frac{1}{2} n^2 \mbox{E}\big( H_n^2\big)\right)\cdot \] A normalidade assintótica de \(\mbox{UCV}\) e \(\mbox{BCV}\) não é óbvia porque o número de termos na soma dupla efetivamente diminui, já que a largura de banda é uma função decrescente de \(n\). Hall usou um argumento da teoria de martingales degenerados para provar esse teorema.

    Scott and Terrell (1987) provaram então que, para uma largura de banda fixa \(h\), as funções \(\mbox{UCV}\) e \(\mbox{BCV}\) eram

    1. ambas assintoticamente normais;

    2. ambas convergiam para \(\mbox{AMISE}(h)\); e

    3. as variâncias assintóticas (verticais) de \(\mbox{UCV}\) e \(\mbox{BCV}\) em \(h\) são \[ \tag{6.73} \dfrac{2 R(\gamma)R(f)}{n^2 h} \qquad \mbox{e} \qquad \dfrac{\sigma_K^8 R(\phi) R(f)}{8n^2h}, \] respectivamente.

    Para kernels da família Beta simétrica da forma (6.27), a variância (vertical) do \(\mbox{UCV}\) é pelo menos 80 vezes maior que a do critério \(\mbox{BCV}\). Essa menor variância vertical sugere que o minimizador real do critério \(\mbox{BCV}\) terá uma variância menor que \(\widehat{h}_{UCV}\).

    As variâncias na equação (6.73) são da ordem de \(O(n^{-9/5})\) se \(h = O(n^{-1/5})\). No entanto, essa rápida taxa de diminuição pode ser explicada pelo fato de que as próprias funções \(\mbox{CV}\) tendem a 0 na taxa \(O(n^{-4/5})\). Assim, a quantidade relevante é o coeficiente de variação (C.V.) \[ C.V.=\dfrac{\sqrt{\mbox{Var}\big(\mbox{UCV}(h) \big)}}{O\big(\mbox{UCV}(h) \big)}=\dfrac{O(n^{-9/10})}{O(n^{-4/5})}=O(n^{-1/10}), \] como afirmado anteriormente.

    O leitor atento notará que o viés não foi contabilizado neste erro no caso do \(\mbox{BCV}\); no entanto, o viés acaba sendo \(O(n^{-1})\) e, portanto, o viés ao quadrado é \(O(n^{-2})\), que é de ordem inferior à da variância.

    Usando um argumento delta descrito abaixo, Scott and Terrell (1987) mostraram que \(\widehat{h}_{UCV}\) e \(\widehat{h}_{BCV}\) convergiam para \(h_{AMISE}\) e eram assintoticamente normais, com variâncias respectivas dadas por \[ \tag{6.74} \dfrac{2 R(f) R\big(\Delta \gamma'(\Delta) \big)}{25 \, n^2 (h^*)^7 \sigma_K^4 R(f'')^2} \qquad \mbox{e} \qquad \dfrac{R(f) R\big(\Delta \phi'(\Delta) \big)}{200\, n^2 (h^*)^7 R(f'')^2} \]

    Observe que se \(h^* = O(n^{-1/5})\), então essas variâncias são \(O(n^{-3/5})\), a partir das quais o coeficiente de variação (horizontal) ainda é de \(O(n^{-1/10})\): \[ C.V.=\dfrac{\sqrt{\mbox{Var}\big(\widehat{h}_{UCV} \big)}}{O\big(\widehat{h}_{UCV} \big)}=\dfrac{O(n^{-3/10})}{O(n^{-1/5})}=O(n^{-1/10})\cdot \] Novamente, para densidades na família Beta simétrica, a variância do \(\mbox{UCV}\) é pelo menos 16 vezes maior que a do \(\mbox{BCV}\). Esses resultados foram confirmados em um estudo de simulação. No entanto, observou-se que o \(\mbox{BCV}\) apresentou desempenho ruim para diversas densidades complexas sem um conjunto de dados muito grande. Essa constatação não é surpreendente, visto que a fórmula do \(\mbox{BCV}\) se baseia no \(\mbox{AMISE}\), enquanto o \(\mbox{MISE}\) exato é a base do \(\mbox{UCV}\).

    É instrutivo descrever a derivação da equação (6.74). Claramente, a largura de banda do \(\mbox{BCV}\) satisfaz \[ \dfrac{\mbox{d}}{\mbox{d}h} \mbox{BCV}(h)\Bigg|_{h=\widehat{h}_{BCV}}=0\cdot \]

    Observando que \[ \Delta'_{ij}=-\dfrac{x_i-x_j}{h^2}=-\dfrac{\Delta_{ij}}{h}, \] a derivada do \(\mbox{BCV}\), conforme definida em (6.71), é igual a \[ -\dfrac{R(K)}{n\, h^2}-\dfrac{\sigma_K^4}{2n^2 h^2}\sum_{i<j} \phi(\Delta_{ij})-\dfrac{\sigma_K^4}{2n^2h}\sum_{i<j} \phi'(\Delta_{ij})\dfrac{\Delta_{ij}}{h}=0 \] ou \[ \sum_{i<j} \Big(\phi(\Delta_{ij})+\Delta_{ij}\phi'(\Delta_{ij}) \Big)\Bigg|_{h=\widehat{h}_{BCV}}=-\dfrac{2nR(K)}{\sigma_K^4}\cdot \]

    Definindo \(\psi(\Delta) = \Delta \phi'(\Delta)\); então, calculando aproximações para os momentos de \(\phi\) e \(\psi\), um cálculo não trivial apresentado na Seção 9 de Scott and Terrell (1987), segue-se que \[ \sum_{i<j} \Big(\phi(\Delta_{ij})+\psi(\Delta_{ij}) \Big)=\mbox{AN}\Big(-2n^2 h^5 R(f''), n^2 h R(f)R(\psi)/2 \Big)\cdot \]

    Portanto, combinando e reorganizando esses dois resultados, \[ -2n^2 R(f'')\widehat{h}^5_{BCV}=\mbox{AN}\Big(-2nR(K)/\sigma_K^4, n^2 h^* R(f)R(\psi)/2 \Big) \] ou \[ \tag{6.75} \widehat{h}^5_{BCV}=\mbox{AN}\left( \dfrac{R(K)}{\sigma_K^4 n R(f'')},\dfrac{h^* R(f)R(\psi)}{8n^2 R(f'')^2}\right)\cdot \]

    Claramente, a média assintótica de \(\widehat{h}^5_{BCV}\) é \((h^*)^5\). A variável aleatória \(\widehat{h}_{BCV}\) é a potência 1/5 daquela dada em (6.75).

    Aplicando o método delta, pode-se concluir que \(\widehat{h}_{BCV}\) é uma \(\mbox{AN}\) com média \(h^*\) e variância que pode ser calculada pela fórmula \[ \mbox{Var}\big(g(h) \big)=\left(\dfrac{\mbox{d}g}{\mbox{d}h} \right)\Bigg|_{h=h^*}^2 \mbox{Var}(h)\cdot \]

    Agora \(g(h)=h^5\) e \(g'(h)=5h^4\), de maneira que \[ \tag{6.76} \mbox{Var}\big(\widehat{h}_{BCV} \big)=\mbox{Var}\big(\widehat{h}_{BCV}^5 \big)/\big( 25(h^*)^8\big)\cdot \] A variância (6.74) surge imediatamente da combinação de (6.75) e (6.76).

    Apesar desses resultados aparentemente favoráveis, o \(\mbox{BCV}\) não se qualifica como um substituto geral para o \(\mbox{UCV}\). O \(\mbox{UCV}\) pode ser mais ruidoso, mas tende a produzir parâmetros de suavização quase não viesados.

    No entanto, há necessidade de um critério auxiliar de validação cruzada, visto que o \(\mbox{UCV}\) é suscetível a certos problemas. Claramente, o \(\mbox{BCV}\) tem suas próprias limitações. Mas, examinando cuidadosamente o trio de parâmetros de suavização sugeridos pelo \(\mbox{UCV}\), \(\mbox{BCV}\) e \(\mbox{OS}\), bem como os formatos das curvas do \(\mbox{UCV}\) e do \(\mbox{BCV}\), boas larguras de banda devem estar disponíveis de forma confiável.

    Na prática, tanto o \(\mbox{UCV}\) quanto o \(\mbox{BCV}\) envolvem computação de ordem \(O(n^2)\) devido às somas duplas. Se as versões gaussianas do kernel forem usadas, o trabalho pode facilmente se tornar proibitivo para \(n > 500\). Esse trabalho pode ser substancialmente reduzido usando uma implementação \(\mbox{ASH}\).

    Os detalhes computacionais são apresentados em Scott and Terrell (1987) e não serão repetidos aqui. Obviamente, existem códigos disponíveis em diversas fontes para realizar esses e outros cálculos. A implementação \(\mbox{ASH}\) é ilustrada na Figura 6.15 para os dados da superfície do aço.

    Figura 6.15: Estimativas de \(\mbox{UCV}\) e \(\mbox{BCV}\) dos dados da superfície do aço \((n = 15.000)\) usando o \(\mbox{ASH}\) triweight. A largura de banda do \(\mbox{UCV}\) foi fixada em 1.2 e 1.3, enquanto a largura de banda do \(\mbox{BCV}\) foi de 1.2 (mostrada): ambas as estimativas foram praticamente idênticas.

    O desvio padrão desses dados é 3.513, de modo que a largura de banda suavizada para o kernel triweight é \(\widehat{h}_{OS} = 1.749\), conforme a equação (6.69). Observe que a estimativa do \(\mbox{UCV}\) é plana em um intervalo relativamente amplo, mesmo com um conjunto de dados tão grande. No entanto, os mínimos dos dois critérios são praticamente idênticos, \(h_{BCV} = 1.2\) e \(h_{UCV} = 1.3\) para o kernel triweight. A estimativa triweight com parâmetro de suavização equivalente \(h = 0.67\) é mostrada na Figura 6.5. Os dados originais foram agrupados em 500 intervalos, de modo que o uso da implementação \(\mbox{ASH}\) das funções de validação cruzada e estimação é natural.


    6.5.1.4 Validação cruzada por Bootstrap


    Taylor (1989) investigou um algoritmo baseado em dados para escolher \(h\) com base em estimativas bootstrap do \(\mbox{MSE}\big(\widehat{f}(x)\big)\) e \(\mbox{MISE}\big(\widehat{f}\big)\).

    A reamostragem bootstrap não é extraída da função de densidade de probabilidade empírica \(f_n\), como no bootstrap ordinário, mas sim a amostra bootstrap \(\{x_1^*, x_2^*,\cdots, x_n^*\}\) é uma amostra aleatória da própria estimativa de densidade kernel candidato \(\widehat{f}(x;h)\).

    Tal reamostragem é chamada de amostra bootstrap suavizada e é discutida mais detalhadamente no Capítulo 9. Considerando \(\mbox{E}_*\) como a esperança em relação à amostra aleatória bootstrap suavizada, Taylor examinou \[ \tag{6.77} \begin{array}{rcl} \mbox{BMSE}_*(x;h) & = & \displaystyle \mbox{E}_*\left( \widehat{f}_*(x;h)-\widehat{f}(x;h)\right)^2 \\[0.8em] & = & \displaystyle \mbox{E}_*\left( \dfrac{1}{n}\sum_{i=1}^n K_h(x-x_i^*) -\dfrac{1}{n}\sum_{i=1}^n K_h(x-x_i)\right)^2 \\[0.8em] \mbox{BMISE}_*(h) & = & \displaystyle \int_x \mbox{BMSE}_*(x;h)\mbox{d}x\cdot \end{array} \]

    O resultado interessante é que, se a reamostragem provém da função de densidade empírica, as quantidades bootstrap na equação (6.77) estimam apenas a variância e não o viés.

    O viés introduzido pelo bootstrap suavizado é precisamente o que é necessário para simular o verdadeiro viés desconhecido para essa escolha do parâmetro de suavização \(h\). Um pouco de variância extra é introduzida pelo bootstrap suavizado, mas é facilmente removida.

    A álgebra envolvida no cálculo do \(\mbox{BMSE}_*\) pode ser desconhecida, mas é simples. Por exemplo, \[ \mbox{E}_*\big(\widehat{f}(x) \big)=\mbox{E}_*\Big( K_h(x-x^*)\Big) \] e \[ \mbox{E}_*\Big( K_h(x-x^*)\Big)=\int K_h(x-y)\widehat{f}(y)\mbox{d}y=\dfrac{1}{n}\sum_{i=1}^n \int K_h(x-y)K_h(y-x_i)\mbox{d}y\cdot \]

    O cálculo de \(\mbox{E}_*\Big( K_h(x-x^*)\Big)^2\) é uma extensão trivial. Para o caso particular do kernel gaussiano, as convoluções indicadas na esperança bootstrap podem ser calculadas em forma fechada.

    Após algumas páginas de trabalho usando o kernel gaussiano e ajustando a variância, Taylor propõe minimizar \[ \tag{6.78} \mbox{BMISE}_*(h)=\dfrac{1}{2nh\sqrt{\pi}}\left(1+\dfrac{\sqrt{2}}{n}\sum_{i<j} \Bigg( \sqrt{2}e^{-\frac{\Delta_{ij}^2}{4}}-\frac{4}{\sqrt{3}}e^{-\frac{\Delta_{ij}^2}{6}} +e^{-\frac{\Delta_{ij}^2}{8}}\Bigg) \right)\cdot \] Taylor demonstra que o \(\mbox{BMISE}_*(h)\) possui a mesma ordem de variância que \(\mbox{UCV}(h)\) e \(\mbox{BCV}(h)\), mas com uma constante menor assintoticamente.

    Na Figura 6.16, as três funções de validação cruzada (CV), com escalas verticais comparáveis, são mostradas graficamente para os dados de queda de neve, juntamente com as estimativas de densidade correspondentes usando o kernel gaussiano. Muitos resultados empíricos comumente observados são representados nesses gráficos.

    Figura 6.16: Algoritmos de validação cruzada com kernel gaussiano e estimativas de densidade para os dados de queda de neve \((n = 63)\). As larguras de banda da validação cruzada são indicadas por setas e a largura de banda suavizada em excesso pela linha tracejada. As estimativas de densidade da validação cruzada não padronizada (\(\mbox{UCV}\)), validação cruzada padronizada (\(\mbox{BCV}\)) e suavizada em excesso são representadas pelas linhas tracejadas, contínuas e pontilhadas, respectivamente.

    As curvas de \(\mbox{BCV}\) e bootstrap são semelhantes, embora \(\mbox{BCV}(h)\) apresente um mínimo mais acentuado. Tanto a função de validação cruzada com viés quanto a função de validação cruzada com bootstrap apresentam mínimos em larguras de banda maiores que a largura de banda suavizada \(h_{OS} = 11.9\). Isso serve para enfatizar a dificuldade em estimar o viés com amostras pequenas. A função de validação cruzada sem viés reflete melhor a dificuldade em estimar o viés com precisão, apresentando uma curva plana em um amplo intervalo próximo ao mínimo. Há alguma evidência visual de três modas nesses dados. Como essas funções CV não levam em conta a natureza de série temporal desses dados e a autocorrelação de primeira ordem é -0.6, uma largura de banda menor é provavelmente justificável Altman (1990).

    Esses cálculos são repetidos para o conjunto de dados do gêiser Old Faithful, que é claramente bimodal, e os resultados são apresentados na Figura 6.17. É interessante notar que tanto a função de validação cruzada com viés quanto a função de validação cruzada com bootstrap apresentam dois mínimos locais.

    Figura 6.17: Algoritmos de validação cruzada com kernel gaussiano e estimativas de densidade para o conjunto de dados de gêiseres \((n = 107)\). As larguras de banda da validação cruzada são indicadas por setas e a largura de banda suavizada por uma linha contínua. As estimativas de densidade da validação cruzada não padronizada (\(\mbox{UCV}\)), validação cruzada padronizada (\(\mbox{BCV}\)) e suavizada são representadas pelas linhas tracejadas, contínuas e pontilhadas, respectivamente.

    Felizmente, apenas um dos mínimos locais é menor que o limite superior de suavização excessiva \(h_{OS} = 0.47\), embora a curva bootstrap mal apresente o segundo mínimo local. A função de validação cruzada com sobresuavização (\(\mbox{UCV}\)) leva a uma largura de banda estreita e a uma estimativa de densidade que parece claramente subsuavizada, dado o tamanho da amostra. Essas observações parecem se repetir com muitos conjuntos de dados “reais” com tamanhos de amostra modestos. A concordância entre um subconjunto desses critérios com comportamentos bastante distintos deve ser levada a sério.


    6.5.1.5 Taxas mais rápidas e validação cruzada PI


    Em uma série de artigos, diversos autores trabalharam para melhorar a convergência relativamente lenta de O(n−1/10 ) dos algoritmos de validação cruzada.


    6.5.1.6 Suavização excessiva restrita


    A suavização excessiva tem sido apresentada estritamente como um meio de delimitar os parâmetros de suavização desconhecidos, escolhendo uma medida de escala, como a amplitude da amostra ou o desvio padrão, para os dados. É fácil encontrar casos em que as larguras de banda suavizadas excessivamente são muito amplas.

    Se a densidade de amostragem for, de fato, a densidade suavizada excessivamente, então a largura de banda ĥOS não é estritamente um limite superior, pois ĥOS varia em torno de h∗. Na maioria dos casos, as duas larguras de banda estarão dentro de alguns pontos percentuais uma da outra.


    6.5.2 Dados multivariados


    A eficácia dos algoritmos de validação cruzada univariada ainda está sob investigação, embora as opções gerais disponíveis sejam razoavelmente bem compreendidas. Os métodos PI podem ser bastante estáveis, enquanto a validação cruzada não enviesada é muito geral e fácil de implementar. Houve progresso nas generalizações multivariadas, embora, novamente, as opções ainda estejam em discussão. Estas são descritas abaixo.


    6.5.2.1 Validação cruzada multivariada


    Em princípio, é simples estender cada um dos algoritmos das seções anteriores para o contexto multivariado. Por exemplo, o algoritmo bootstrap é facilmente estendido, assim como as expressões analíticas para \(\mbox{BCV}\) e \(\mbox{UCV}\) (Sain, Baggerly, and Scott 1994). Exemplos de \(\mbox{UCV}\), \(\mbox{BCV}\) e \(\mbox{BMISE}^* (h_1,h_2)\) baseados em 1500 pontos \(N(\pmb{0}_2,\pmb{I}_2)\) são mostrados na Figura 6.24. Todos são razoavelmente próximos neste caso.

    Figura 6.24: \(\mbox{MISE}(h_x, h_y)\) estimado usando \(\mbox{UCV}\), \(\mbox{BCV}\) e os algoritmos de bootstrap em 1500 pontos \(N(\pmb{0}_2,\pmb{I}_2)\). Os dois pontos em cada diagonal representam o \(h^*\) e as larguras de banda suavizadas. O ponto que localiza o minimizador de cada critério está abaixo da diagonal.

    A Figura 6.25 mostra os mesmos critérios aplicados ao conjunto de dados log-lipídicos padronizado. No entanto, a estimativa \(\mbox{BMISE}^∗(h_1,h_2)\) para os valores lipídicos \((n = 320)\) tem seu mínimo além da regra de referência normal e das larguras de banda excessivamente suavizadas. As larguras de banda de \(\mbox{UCV}\) e \(\mbox{BCV}\) são semelhantes, embora a \(\mbox{BCV}\) seja (surpreendentemente) um pouco mais agressiva. Ambos mostram o segundo e possivelmente o terceiro moda/pico, enquanto a estimativa Bootstrap é unimodal.

    Figura 6.25: O mesmo critério da Figura 6.24 foi aplicado ao conjunto de dados do logaritmo do lipídio padronizado \((n = 320)\), juntamente com as estimativas de kernel correspondentes.

    Em geral, pode-se esperar que o desempenho de todos os algoritmos de validação cruzada multivariada seja mais assintótico e, portanto, mais lento para convergir na prática. Manter um parâmetro de suavização para cada direção de coordenada é importante. Mesmo que os dados tenham sido padronizados, é improvável que \(h_i = h_j\) seja satisfatório Wand and Jones (1991a)

    Menos claro é se o ganho obtido com o uso de uma matriz de suavização completa, que permite contornos elípticos, é suficiente para compensar o aumento drástico no número de parâmetros que devem ser estimados por validação cruzada. Em muitos ou na maioria dos casos, os ganhos não são obtidos na prática. Além disso, uma covariância completa pode degenerar em uma forma singular, que pode ter probabilidade “infinita”. Valores iniciais diferentes podem evitar esse problema.

    A estimativa de densidade também é bastante difícil se a nuvem de dados for altamente assimétrica ou se concentrar (quase) em uma variedade de dimensão inferior. Este tópico é abordado em detalhes no Capítulo 7. No entanto, transformações marginais usando a escada de Tukey, por exemplo, são sempre recomendadas. Mas, assim como no cenário univariado, uma escolha absolutamente ótima da largura de banda não é crítica para fins exploratórios.


    6.5.2.2 Larguras de banda de sobresuavização multivariada


    A extensão da largura de banda de sobresuavização para d foi resolvida por Terrell (1990). A densidade multivariada mais simples é esfericamente simétrica. Assim, a formulação geral do kernel é necessária com a restrição de que todas as larguras de banda marginais sejam iguais. Por simetria, encontrar a forma da densidade sobresuavizada multivariada ao longo do eixo x, por exemplo, será suficiente.


    6.6 Suavização adaptativa


    6.6.1 Introdução ao kernel variável


    Considere o estimador kernel fixo multivariado \[ \widehat{f}(\pmb{x})=\dfrac{1}{nh^d}\sum_{i=1}^n K\Bigg( \dfrac{\pmb{x}-\pmb{x}_i}{h}\Bigg)=\dfrac{1}{n}\sum_{i=1}^n K_h(\pmb{x}-\pmb{x}_i)\cdot \]

    O estimador adaptativo mais geral dentro dessa estrutura simples permite que a largura de banda \(h\) varie não apenas com o ponto de estimação, mas também com a realização particular da densidade desconhecida \(f\): \[ h \leftarrow h(\pmb{x},\pmb{x}_i,\{\pmb{x}_j\}) \approx h(\pmb{x},\pmb{x}_i,f)\cdot \]

    A segunda forma indica que, assintoticamente, a porção da fórmula da largura de banda adaptativa que depende da amostra completa pode ser representada como uma função da densidade verdadeira. Além disso, pode-se assumir que a função de largura de banda adaptativa ótima, \(h(\pmb{x},\pmb{x}_i)\), é suave e varia lentamente.

    Assim, para amostras finitas, será suficiente considerar duas abordagens distintas para a estimação adaptativa da função de densidade: \[ \tag{6.83} \widehat{f}_1(\pmb{x})=\dfrac{1}{n}\sum_{i=1}^n K_{h_{\pmb{x}}}(\pmb{x}-\pmb{x}_i), \qquad h_{\pmb{x}}=h(\pmb{x},\pmb{x}_i,f) \] ou \[ \tag{6.84} \widehat{f}_2(\pmb{x})=\dfrac{1}{n}\sum_{i=1}^n K_{h_i}(\pmb{x}-\pmb{x}_i), \qquad h_i=h(\pmb{x}_i,\pmb{x}_i,f)\cdot \]

    No primeiro caso, uma largura de banda fixa é usada para todos os \(n\) pontos de dados, mas essa largura de banda fixa muda para cada ponto de estimativa \(\pmb{x}\). No segundo caso, uma largura de banda diferente é escolhida para cada \(\pmb{x}_i\) e, em seguida, aplicada para estimar a densidade globalmente.

    Cada uma é justificada assintoticamente pela suposição de suavidade local em \(h(\pmb{x},\pmb{x}_i,f)\), uma vez que, assintoticamente, apenas os pontos de dados em uma pequena vizinhança de \(\pmb{x}\) contribuem para o valor da densidade ali.

    Presumivelmente, todas as larguras de banda ótimas nessa vizinhança são suficientemente próximas para que o uso de apenas um valor seja adequado. A escolha de \(\widehat{f}_1\) ou \(\widehat{f}_2\) depende das dificuldades práticas em especificar a função de largura de banda adaptativa. Para amostras pequenas, pode-se esperar alguma diferença de desempenho entre os dois estimadores. Jones (1990) apresentou uma demonstração gráfica útil das diferenças entre os dois estimadores adaptativos.

    Exemplos de \(\widehat{f}_1\) incluem o estimadores \(k\)-NN de Loftsgaarden and Quesenberry (1965) (ver Seção 6.4.1) com \(h_{\pmb{x}}\) igual à distância até o \(k\)-ésimo ponto de amostragem mais próximo: \[ \tag{6.85} h_{\pmb{x}}=d_k(\pmb{x},\{\pmb{x}_i\})\approx \left(\dfrac{k}{n \, V_df(\pmb{x})} \right)^{1/d}, \] onde a distância estocástica é substituída por uma fórmula simples semelhante a um histograma.

    A segunda forma foi introduzida por Breiman, Meisel, and Purcell (1977), que sugeriram a escolha de \[ \tag{6.86} h_i = h\times d_k\big(\pmb{x}_i,\{\pmb{x}_1,\pmb{x}_2,\cdots,\pmb{x}_n\} \big)\approx h\times \left( \dfrac{k}{n V_d f(\pmb{x})}\right)^{1/d}\cdot \] A semelhança entre essas duas propostas específicas é evidente. No entanto, o foco no uso da distância \(k\)-NN é simplesmente uma questão de conveniência.

    Quando estimado de forma ótima ponto a ponto, o primeiro método fornece a melhor estimativa possível assintoticamente, pelo menos do ponto de vista do \(\mbox{MISE}\). Contudo, o estimador não é, por construção, uma função de densidade.

    Por exemplo, é fácil ver que o estimador \(k\)-NN integra para \(\infty\) (veja o Exercício 34). O segundo estimador, por outro lado, é, por construção, um estimador de densidade genuíno para kernels não negativos.

    Na prática, estimadores adaptativos nas formas muito gerais dadas nas equações (6.83) e (6.84) podem ser muito difíceis de especificar. A “função de suavização adaptativa” \(h_{\pmb{x}}\) para \(\widehat{f}_1\) é \(\infty\)-dimensional. A especificação para \(\widehat{f}_2\) é um pouco mais fácil, já que o “vetor de suavização adaptativa” \(\{h_i\}\) é apenas \(n\)-dimensional. Como de costume, as escolhas “corretas” dessas quantidades dependem de um conhecimento mais aprofundado das derivadas desconhecidas da função de densidade.

    Na prática, os estimadores adaptativos tornam-se viáveis reduzindo-se significativamente a dimensão da função de suavização adaptativa. Um exemplo simples é incorporar a função de distância \(d_k(\cdot,\cdot)\) como nas equações (6.85) e (6.86). Abramson (1982) propôs uma variação da fórmula de Breiman et al. (6.86): \[ h_i=h/\sqrt{f(\pmb{x_i})}, \qquad \forall \pmb{x}_i \in \mathbb{R}^d\cdot \] Em um artigo complementar, Abramson (1982b) demonstra que o uso de uma estimativa piloto não adaptativa para f é adequado. Observe que essas duas propostas concordam quando d = 2, onde Breiman et al. (1977) descobriram empiricamente que sua fórmula funcionava bem.

    No contexto univariado, a ideia simples de aplicar um estimador kernel fixo aos dados transformados e, em seguida, realizar a transformação inversa também se enquadra na segunda categoria. Se \(u\) for uma função monótona suave selecionada de uma família de transformações como a de Box and Cox (1964), então, quando a estimativa kernel fixo de \(\omega = u(x)\) for retransformada para a escala original, o efeito é especificar implicitamente o valor de \(h_i\), pelo menos assintoticamente.

    A abordagem de transformação demonstrou capacidade de lidar bem com dados assimétricos e, em menor grau, com dados curtos simétricos (ver Wand and Jones 1991b). A técnica de transformação não funciona tão bem com dados multimodais.

    Em cada caso, a potencial instabilidade do estimador adaptativo foi significativamente reduzida pelo “enrijecimento” da função ou vetor de suavização. Isso pode ser visto explicitamente contando o número de parâmetros de suavização \(s\) que devem ser especificados:

    1. \(s = 1\) para \(k\)-NN \((k)\);

    2. \(s = 2\) para Breiman et al. \((h,k)\);

    3. \(s = 2\) para Abramson \((h,h_{pilot})\); e

    4. \(s = 2\), 3 ou 4 para a abordagem de transformação, \(h\) mais 1-3 parâmetros que definem a família de transformação.

    Os aspectos teóricos e práticos do problema adaptativo são investigados mais a fundo. Sabe-se muito mais sobre os primeiros do que sobre os últimos.


    6.6.2 Suavização adaptativa univariada


    6.6.2.1 Limites de melhoria


    Considere um estimador adaptativo pontual com um kernel de ordem \(p\): \[ \widehat{f}_1(x)=\dfrac{1}{n\, h(x)}\sum_{i=1}^n K\Bigg(\dfrac{x-x_i}{h(x)} \Bigg), \qquad \mbox{permitindo} \qquad h(x)=h_{x}\cdot \]

    As propriedades \(\mbox{AMSE}\) pontuais deste estimador podem ser reconhecidas no Teorema 6.3 para um kernel de ordem \(p\): \[ \mbox{AV}(x)=\dfrac{R(K) f(x)}{n\, h(x)} \qquad \mbox{e} \qquad \mbox{ASB}(x)=\dfrac{\mu_p^2 f^{(p)}(x)^2}{(p!)^2}h(x)^{2p}, \] a partir do qual pode-se obter o \(\mbox{AMSE}(x)\) pontual ótimo: \[ h^*= \left( \dfrac{(p!)^2 R(K)f(x)}{2p\,\mu_p^2\, f^{(p)}(x)^2 }\right) n^{-1/(2p+1)} \] e, portanto, \[ \tag{6.87} \mbox{AMSE}^*(x)=\dfrac{2p+1}{2p}\left( \dfrac{2p \, \mu_p^2 \, R(K)^{2p} f(x)^{2p} f^{(p)}(x)^2 }{(p!)^2}\right)^{1/(2p+1)} n^{-2p/(2p+1)}\cdot \]

    Lembre-se de que o \(\mbox{MISE}\) acumula erros pontuais. Assim, acumular os erros pontuais mínimos obtidos usando \(h^*(x)\) fornece o limite inferior assintótico para o \(\mbox{AMISE}\) adaptativo: \[ \tag{6.88} \begin{array}{rcl} \mbox{AAMISE}^* & = & \displaystyle \int_{-\infty}^\infty \mbox{AMSE}^*(x)\mbox{d}x \\[0.8em] & = & \displaystyle \dfrac{2p+1}{2p}\left( \dfrac{2p \, \mu_p^2 \, R(K)^{2p}}{(p!)^2}\right)^{1/(2p+1)}\times n^{-2p/(2p+1)} \int_{-\infty}^\infty \Big( f(x)^{2p} f^{(p)}(x)^2\Big)^{1/(2p+1)} \mbox{d}x\cdot \end{array} \]

    Comparando as equações (6.23) e (6.88), conclui-se que o limite para a melhoria de um estimador kernel adaptativo é \[ \tag{6.89} \dfrac{\mbox{AAMISE}^*}{\mbox{AMISE}^*} = \dfrac{\displaystyle \int_{-\infty}^\infty \Big( f(x)^{2p} f^{(p)}(x)^2\Big)^{1/(2p+1)}\mbox{d}x}{\displaystyle \int_{-\infty}^\infty \Big( f^{(p)}(x)^2\Big)^{1/(2p+1)} \mbox{d}x}\cdot \]

    Uma aplicação da desigualdade de Jensen à quantidade \[ \mbox{E}\left( f^{(p)}(X)^2 \Big/ f(X)\right)^{1/(2p+1)} \] mostra que a razão em (6.89) é sempre \(\leq 1\) (ver Exercício 35). Na Tabela 6.4, essa razão de limite inferior é calculada numericamente para as densidades normal e de Cauchy.

    Tabela 6.4: Razão entre \(\mbox{AAMISE}^*\) e \(\mbox{AMISE}^*\) para duas densidades comuns em função da ordem do kernel.

    Observe que o potencial de adaptabilidade diminui para kernels de ordem superior se os dados forem normais, mas o oposto ocorre para dados de Cauchy. A tabela fornece mais evidências da relativa facilidade na estimativa da densidade normal. Rosenblatt derivou (6.88) no caso de kernel positivo \(p = 2\).

    No caso de um kernel positivo \(p = 2\), a malha adaptativa assintoticamente ótima é, na verdade, igualmente espaçada quando \[ \tag{6.90} \dfrac{f''(x)^2}{f(x)}=c \qquad \mbox{o qual implica} \qquad f(x)=\dfrac{c}{144}(x-a)^4, \] onde \(a\) é uma constante arbitrária (ver Seção 3.2.8.3).

    Assim, o espaço nulo para estimadores kernel ocorre em intervalos onde a densidade é uma função quártica pura. A junção de segmentos puros da forma (6.90), garantindo que \(f\) e \(f'\) sejam contínuos, implica que \(f\) é monótona; portanto, não existe uma densidade adaptativa nula completa em \(C^1\) ou \(C^2\), a menos que kernels de fronteira sejam introduzidos.


    6.6.2.2 Estimadores de vizinhos mais próximos


    Usando o valor assintótico para o parâmetro de suavização adaptativa da equação (6.85) com um kernel positivo \((p = 2)\) e \(d = 1\), \[ \tag{6.91} h(x)=\dfrac{k}{2n \, f(x)}, \] o viés quadrático integrado assintótico adaptativo é dado por \[ \mbox{AAISB}(k)=\int_x \mbox{AISE}(x)\mbox{d}x =\int_x \dfrac{1}{4}h(x)^4 f''(x)^2 \mbox{d}x=\dfrac{k^4}{64 \, n^4} \int_x \dfrac{f''(x)^2}{f(x)^4}\mbox{d}x\cdot \]

    Surpreendentemente, observa-se facilmente que a última integral diverge para densidades simples como a normal. Essa divergência não implica que o viés seja infinito para amostras finitas, mas indica que nenhuma escolha de \(k\) pode proporcionar uma relação satisfatória entre variância e viés.

    A explicação é bastante simples: assintoticamente, o suavização adaptativa ótima depende não apenas do nível de densidade, como na equação (6.91), mas também da curvatura, como na equação (6.87). Assim, a regra simples dada pela equação (6.91) não representa uma melhoria em relação à estimação não adaptativa.

    Esse fenômeno de pior desempenho será observado repetidamente com muitos procedimentos adaptativos ad hoc simples, pelo menos assintoticamente. Alguns proporcionam ganho significativo para amostras pequenas ou certas densidades.

    Outras abordagens, como a transformação, incluem a configuração de largura de banda fixa como um caso especial e, portanto, não precisam apresentar desempenho significativamente pior assintoticamente.


    6.6.2.3 Estimadores adaptativos por ponto de amostragem


    Considere a segunda forma para um estimador adaptativo com diferentes parâmetros de suavização nos pontos de amostragem: \[ \tag{6.92} \widehat{f}_2(x)=\dfrac{1}{n}\sum_{i=1}^n \dfrac{1}{h(x_i)} K\Bigg( \dfrac{x-x_i}{h(x_i)}\Bigg), \qquad \mbox{permitindo} \qquad h(x_i)=h_i\cdot \]

    G. R. Terrell and Scott (1992a) provaram, sob certas condições, que \[ \mbox{AV}\big(\widehat{f}_2(x)\big)=\dfrac{f(x)R(K)}{n\, h(x)} \] e \[ \mbox{ASB}\big(\widehat{f}_2(x)\big)=\left( \dfrac{\big(h(x)^p f(x) \big)^{(p)}}{p!}\right)^2, \] onde o sobrescrito \((p)\) indica uma derivada de ordem \(p\).

    Abramson (1982) propôs o que é óbvio a partir desta expressão para o viés, ou seja, que quando \(p = 2\), a escolha \[ \tag{6.93} h(x)=h\Big/\sqrt{f(x)} \] implica que o termo de viés de segunda ordem em \(\mbox{ASB}\) desaparece! O viés é, na verdade, \[ \dfrac{1}{4!}\big(h(x)^4 f(x) \big)^{(iv)}=\dfrac{1}{24}h^4 \left(\dfrac{1}{f(x)} \right)^{(iv)}, \] como demonstrado por Silverman (1986).

    Esse viés de quarta ordem geralmente é reservado para kernels negativos de \(p=4\) ordem e aparentemente contradiz o resultado clássico de Farrell (1972) sobre as melhores taxas de viés com kernels positivos.

    De fato, G. R. Terrell and Scott (1992a) forneceram um exemplo que ilustra o comportamento real com dados normais \((f = \phi)\): \[ \widehat{f}_2(x)=\dfrac{1}{n}\sum_{i=1}^n \dfrac{\sqrt{\phi(x_i)}}{h} K\Bigg(\dfrac{(x-x_i)\sqrt{\phi(x_i)}}{h} \Bigg), \] do que se conclui que \[ \tag{6.94} \mbox{E}\big(\widehat{f}_2(x) \big) = \int \dfrac{\sqrt{\phi(t)}}{h}K\Bigg(\dfrac{(x-t)\sqrt{\phi(t)}}{h} \Bigg),\phi(t)\mbox{d}t, \] com uma expressão similar para a variância.

    O \(\mbox{MISE}\) adaptativo exato é difícil de obter, mas pode ser calculado numericamente para escolhas específicas de \(h\) e \(n\). Os autores mostraram que o \(\mbox{MISE}\) adaptativo exato era metade do \(\mbox{MISE}\) do melhor estimador de largura de banda fixa quando \(n< 200\); no entanto, o estimador de largura de banda fixa foi superior para \(n > 20.000\).

    Essa descoberta sugere que o procedimento não possui \(\mbox{MISE}\) de ordem \(O(n^{-8/9})\). De fato, pode-se demonstrar a origem da diferença com Silverman (1986) e Hall and Marron (1988). Considere uma fórmula assintótica para o \(\mbox{MSE}\) do estimador em \(x = 0\), sem perda de generalidade. O ponto sutil é evidenciado mais claramente pela escolha do kernel boxcar \(K = U(-1, 1)\) de modo que a equação (6.94) se torne \[ \tag{6.95} \mbox{E}\big(\widehat{f}(0) \big) =\dfrac{1}{2h}\int \phi(t)^{3/2} \left( \pmb{I}_{[-1,1]} \Bigg(\dfrac{t\sqrt{\phi(t)}}{h}\Bigg)\right)\mbox{d}t\cdot \]

    Normalmente, os limites de integração se estenderiam de \(-h\) a \(h\). No entanto, uma análise mais detalhada mostra que o argumento do kernel não é monotonicamente crescente e, na verdade, tende a zero quando \(|t|\to\infty\).

    Assim, a integral na equação (6.95) abrange três intervalos, que chamaremos de \((-\infty,-b)\), \((-a,a)\) e \((b,+\infty)\), onde \(a\) e \(b\) são soluções da equação \[ \tag{6.96} \dfrac{t\sqrt{\phi(t)}}{h}=1\cdot \]

    Defina \(c = (2\pi)^{1/4} h\). Então a equação (6.96) assume duas formas que fornecem aproximações suficientes para os extremos do intervalo \(a\) e \(b\): \[ t \exp\big(-t^2/4 \big)=c \qquad \mbox{implica} \qquad a\approx c+ c^3/4+5c^5/32,\\ \log(t)-t^2/4=\log(c) \qquad \mbox{implica} \qquad b\approx \Big( -4\log(c)+4\log\big(\sqrt{-4\log(c)}\big)\Big)^{1/2}\cdot \]

    Agora, tomando a série de Taylor do integrando na equação (6.95) e integrando em \((-a,a)\), obtemos: \[ \mbox{E}\big(\widehat{f}(0) \big)=\phi(0)+(2\pi)^{1/2} h^4/40+O(h^6), \] o que fornece o viés previsto de \(O(h^4)\).

    No entanto, a contribuição para o viés dos dois intervalos restantes totaliza \[ 2\times \dfrac{1}{2h}\int_{-\infty}^{-b} \phi(x)^{3/2}\mbox{d}x =\left(\dfrac{2}{9\pi}\right)^{1/4} \dfrac{1}{h}\Phi\Big(-b\sqrt{3/2} \Big), \] que, usando a aproximação para \(b\) e a aproximação da cauda \(\Phi(x)\approx -\phi(x)/x\) para \(x<<0\), iguala \[ h^2\Big/ \Big(24\big(\log( (2\pi)^{1/4} h) \big)^2 \Big)\cdot \]

    Assim, as caudas exercem uma influência indevida na estimativa no meio e destroem o ganho aparente em viés. Com um kernel mais suave, o mesmo efeito é observado, mas não tão claramente. Abramson reconheceu esse problema prático e sugeriu impor um limite superior a \(h_i\) “cortando” o estimador piloto na equação (6.93) longe de zero.

    Outros autores não consideraram essa sugestão. No entanto, a ineficiência assintótica não anula as boas propriedades em pequenas amostras observadas por Abramson (1982), Silverman (1986) e Worton (1989). Essa mesma análise pode ser aplicada à proposta original de Breiman, Meisel, and Purcell (1977). A contribuição das caudas para o viés acaba sendo da ordem \(O\big(h/\log(h)\big)\).

    Esses autores observaram que, apesar do excelente desempenho bivariado empírico, o desempenho univariado era ruim. Essa baixa taxa de viés ajuda a explicar essa observação.


    6.6.2.4 Aprimoramento de dados


    O viés de uma estimativa kernel positiva \(\widehat{f}(x)\), é controlado pela segunda derivada nesse ponto (ver equação (6.16)). Assim, as modas geralmente serão subestimados e as antimodas sempre serão superestimados. Mesmo se forem empregados kernels negativos de ordem superior, esse mesmo fenômeno será observado. Lembre-se de que o método do gradiente minimiza uma função movendo-se na direção do gradiente negativo.

    Samiuddin and El-Sayyad (1990) propuseram ajustar os dados em direção ao pico mais próximo, seguindo o gradiente positivo. Especificamente, no cenário univariado em que \(X_i\sim f(x)\), eles propõem usar uma estimativa kernel nos dados ajustados \[ \tag{6.97} \widetilde{x}_i = x_i+\dfrac{h^2}{2}\dfrac{f'(x_i)}{f(x_i)}, \qquad \mbox{para} \qquad i=1,\dots,n \] que, no caso gaussiano seria \[ \tag{6.98} \widetilde{x}_i=x_i-\dfrac{h^2}{2}\dfrac{x_i-\mu}{\sigma^2}, \qquad \mbox{para} \qquad x\sim N(\mu,\sigma^2), \] onde assumimos que o kernel é escalado de forma que \(\sigma_K = 1\).

    Note que, para uma densidade normal, a equação (6.98) pode ser reescrita como \((\widetilde{x}_i-\mu) = (x_i-\mu)\big(1-h^2/(2\sigma^2)\big)\) ; portanto, os dados são comprimidos em direção à média (e à moda) pelo fator \(\frac{1}{2}\big(2-h^2/\sigma^2\big)\). Versões amostrais da equação (6.97) e também de forma mais geral foram consideradas por Choi and Hall (1999) e Hall and Minnotte (2002).

    Para investigar a eficácia desta proposta, focamos na eeperança do estimador kernel usando os dados modificados: \[ \tag{6.99} \widehat{f}(x)=\dfrac{1}{n}\sum_{i=1}^n K_h(x-\widetilde{x}_i)=\dfrac{1}{n}\sum_{i=1}^n \dfrac{1}{h}K\left(\dfrac{x-x_i-\frac{1}{2}h^2 f'(x_i)/f(x_i)}{h} \right), \] onde a densidade verdadeira é usada no ajuste (6.97).

    Suponha que o kernel \(K\) seja simétrico, \(K(-\omega) = K(\omega)\), e que \(\sigma_K = 1\). Usando a mudança de variáveis \(\omega = (x-y)/h\) e uma Série de Taylor no kernel, a esperança é dada por \[ \begin{array}{rcl} \mbox{E}\big(\widehat{f}(x)\big) & = & \displaystyle \int \dfrac{1}{h}K\left(\dfrac{x-y-\frac{1}{2}h^2 f'(y)/f(y)}{h} \right)f(y)\mbox{d}y\\[0.8em] & = & \displaystyle \int K\left(\omega-\dfrac{h}{2} \dfrac{f'(x-h\omega)}{f(x-h\omega)} \right)f(x-h\omega)\mbox{d}y \\[0.8em] & = & \displaystyle \int f(x-h\omega)\left( K(\omega)-\dfrac{h}{2}\dfrac{f'(x-h\omega)}{f(x-h\omega)}K'(\omega)+\dfrac{h^2}{8}\dfrac{f'(x-h\omega)^2}{f(x-h\omega)^2}K''(\omega) \right)\mbox{d}\omega \\[0.8em] & = & \displaystyle \int K(\omega)f(x-h\omega)\mbox{d}\omega -\dfrac{h}{2}\int K'(\omega)f'(x-h\omega)\mbox{d}\omega+\dfrac{h^2}{8}\int K''(\omega)\dfrac{f'(x-h\omega)^2}{f(x-h\omega)^2}\mbox{d}\omega \\[0.8em] & = & \displaystyle \left(f(x)+\dfrac{1}{2}h^2 f''(x)+O(h^4) \right) - \dfrac{h}{2}\left(h f''(x)+O(h^3) \right) + \dfrac{h^2}{8}O(h^2) \\[0.8em] & = & \displaystyle f(x)+\dfrac{1}{2}h^2 f''(x)-\dfrac{1}{2}h^2 f''(x)+O(h^4) = f(x)+O(h^4)\cdot \end{array} \]

    O resultado é que o viés da proposta de Samiuddin e El-Sayyad é de ordem superior \(O(h^4)\), mesmo com um kernel positivo. Na antepenúltima expressão, a primeira integral é a expansão de viés já conhecida.

    A segunda integral pode ser aproximada aplicando-se uma Série de Taylor a \(f(x-h\omega)\): \[ \int K'(\omega) \left( \sum_{\ell=0} (-h\omega)^\ell f^{(1+\ell)}(x)/\ell! \right)\mbox{d}\omega = 0-h f''(x)+0+O(h^3), \] visto que \(\displaystyle \int K'(\omega)\mbox{d}\omega = 0\), \(\displaystyle \int \omega K'(\omega)\mbox{d}\omega = -1\), \(\displaystyle \int \omega^2 K'(\omega)\mbox{d}\omega = 0\) e \(\displaystyle \int w^3 K'(\omega)\mbox{d}\omega = -3\sigma_K^2\) .

    Pouco se ganha ao derivar a expressão explícita do viés, que é bastante complexa (ver Hall and Minnotte 2002). Os termos de variância principais acabam sendo os mesmos que para o estimador kernel ordinário. Embora essa derivação ajuste os dados usando \(f(x)\), que é desconhecida, o algoritmo amostral explicado por Hall and Minnotte (2002) atinge qualquer ordem superior desejada.

    Como já vimos diversas vezes, a busca por um algoritmo que reduza o viés por meio de adaptação ou ajuste local resulta, na verdade, em um viés de ordem superior. Hall and Minnotte (2002) e outros autores relatam exemplos promissores com dados reais e algumas simulações.

    No entanto, a aplicação global do ajuste de dados (6.97) apresenta uma deficiência evidente, como pode ser observado ao analisar sua ação sobre a densidade da mistura na equação (3.78). A partir da equação (6.98), para dados normais padrão, o ajuste é \(-h^2 x_i/2\), enquanto para dados \(N(3, 1/3^2)\), o ajuste é \(-9h^2 (x_i-3)/2\), que é, na verdade, 27 vezes maior em unidades padrão.

    Intuitivamente, o ajuste para os dados de mistura mais estreitos deveria ter um terço da largura, visto que o desvio padrão é 1/3. Uma possível solução seria incorporar a magnitude da segunda derivada na moda local. Ao examinar os estudos de caso publicados, muitos apresentavam a mesma curvatura nas modas e o efeito passou despercebido. Contudo, essa ideia se mostra promissora com uma implementação de estimação em múltiplos estágios para encontrar a localização das modas e antimodas, bem como a escala local.

    Também pode ser aconselhável usar parâmetros de suavização diferentes nas equações (6.97) e (6.99) e determinar ambos os parâmetros de suavização por meio de validação cruzada não enviesada. Em geral, cada região modal pode exigir seu próprio parâmetro de suavização.


    6.6.3 Procedimentos adaptativos multivariados


    Os procedimentos adaptativos multivariados contêm algumas características interessantes e únicas. Esses resultados estão disponíveis em mais detalhes em G. R. Terrell and Scott (1992b).


    6.6.3.1 Adaptação pontual


    Seja \(\nabla^2 f(\pmb{x})\), que é a matriz das segundas derivadas parciais de \(f\) em \(\pmb{x}\), seja denotada por \(S_{\pmb{x}}\). Da equação (6.48), o viés assintótico pontual é \[ \tag{6.100} \mbox{AB}(\pmb{x})=\dfrac{1}{2}h^2 \mbox{tr}\big(A^\top S_{\pmb{x}} A \big)=\dfrac{1}{2}h^2 \mbox{tr}\big( AA^\top S_{\pmb{x}}\big)\cdot \]

    No contexto univariado, o viés é controlado inteiramente pela escala do kernel, enquanto no contexto multivariado, a forma do kernel também pode ser usada para controlar o viés.

    A importância da forma na minimização de \(\mbox{tr}\big(A^\top S_{\pmb{x}} A\big)\) depende das propriedades da matriz \(S_{\pmb{x}}\).

    Caso 1 \(S_{\pmb{x}}\) é definida positiva ou negativa. Como \(H\) e, portanto, \(A\) tem posto completo por hipótese, então \(A^\top S_{\pmb{x}} A\) também é definida positiva ou negativa. Como a soma e o produto dos autovalores de uma matriz definida são iguais ao seu traço e determinante, respectivamente, a matriz com traço e determinante mínimos (absolutos) iguais a 1 tem todos os seus autovalores iguais.

    Portanto, a matriz \(A\) deve ser escolhida para satisfazer \[ A A^\top S_{\pmb{x}}=|S_{\pmb{x}}|^{1/d} \pmb{I}_d\cdot \]

    Observe que, com essa escolha, a matriz do lado direito tem o mesmo determinante que \(S_{\pmb{x}}\), e todos os autovalores de \(AA^\top S_{\pmb{x}}\) são iguais a \(|S_{\pmb{x}}|^{1/d}\). Portanto, o melhor \[ \mbox{tr}\big(A^\top S_{\pmb{x}} A\big) = d|S_{\pmb{x}}|^{1/d}\cdot \] O \(\mbox{MSE}\) assintótico pontual de \(\widehat{f}(\pmb{x})\) segue das equações (6.49) e (6.100) e pode ser otimizado para produzir \[ h^*= \left( \dfrac{f(\pmb{x})R(K)}{n \,d \,|S_{\pmb{x}}|^{2/d}}\right)^{1/(d+4)} \] de maneira que \[ \mbox{AMSE}^*(\pmb{x})= \left( \dfrac{d+4}{4d}\right)^{2(d+2)/(d+4)} \Big(f(\pmb{x})R(K)\sqrt{|S_{\pmb{x}}|}\Big)^{4/(d+4)} n^{-4/(d+4)}\cdot \]

    Caso II A densidade tem formato de sela em \(\pmb{x}\); isto é, a matriz \(S_{\pmb{x}}\) possui autovalores positivos e negativos. A densidade é curvada para cima em algumas direções e para baixo em outras.

    Nesse caso, é possível construir a matriz \(A\) de modo que a soma dos autovalores de \(AA^\top S_{\pmb{x}}\) seja igual a 0 (ver G. R. Terrell and Scott 1992b). Assim, os termos de viés de ordem \(h^2\) desaparecem, com os termos de ordem superior dominando.

    Portanto, regiões onde a densidade tem formato de sela não contribuem assintoticamente para o \(\mbox{AAMISE}\) em comparação com regiões onde a densidade é definida.

    Quão comuns são as regiões com formato de sela? Considere a densidade normal multivariada \(N(\pmb{0}_d,\pmb{I}_d)\). Então, o gradiente e a matriz Hessiana de \(f(\pmb{x})\) são \[ \nabla f(\pmb{x})=-\nabla f(\pmb{x}) \qquad \mbox{e} \qquad \nabla^2 f(\pmb{x})=S_{\pmb{x}} =f(\pmb{x})\big(\pmb{x}\pmb{x}^\top-\pmb{I}_d \big)\cdot \]

    Existem \(d-1\) autovalores da matriz \(\big(\pmb{x}\pmb{x}^\top-\pmb{I}_d \big)\) iguais a -1 e um igual a \(\pmb{x}\pmb{x}^\top-1\). Portanto, a densidade normal multivariada é negativa definida quando \(||\pmb{x}|| < 1\) e em forma de sela em todos os pontos fora da esfera unitária.

    Dado que a fração da massa de probabilidade dentro da esfera unitária diminui à medida que a dimensão aumenta, o potencial significado prático dessa descoberta parece promissor. Finalmente, observe que exatamente na esfera unitária onde \(||\pmb{x}|| = 1\), um dos autovalores se anula. Isso leva ao caso final.

    Caso III \(S_{\pmb{x}}\) é semidefinido com pelo menos um autovalor zero. A densidade é plana em certas direções. Não há contribuição para o viés nessas direções e, portanto, a contribuição do viés é novamente de ordem inferior à do Caso I. O problema, nesse ponto, pode ser reduzido a uma dimensão inferior por projeção, se desejado.


    6.6.3.2 Adaptação global


    O problema de selecionar o melhor estimador kernel fixo adaptativo global foi formulado por Deheuvels (1977a), Deheuvels (1977b), que caracterizou a solução na forma de equação diferencial.

    O critério global a ser otimizado é \[ \mbox{AAMISE}(h,A)=\dfrac{R(K)}{n\, h^d}+\dfrac{1}{4}h^4 \int_{\mathbb{R}^d} \mbox{tr}\big(A^\top S_{\pmb{x}}A\big) \mbox{d}\pmb{x}\cdot \]

    Minimizar em relação a \(h\) e \(A\) separadamente resulta em \[ \mbox{AAMISE}^* = \left(\min_A \int_{\mathbb{R}^d} \mbox{tr}^2\big(A^\top S_{\pmb{x}}A\big) \mbox{d}\pmb{x} \right)^{d/(d+4)} \left(\dfrac{(d+4)R(K)}{4nd} \right)^{4/(d+4)}\cdot \]

    Para dados \(N(\pmb{0}_d,\pmb{I}_d)\), a vantagem do esquema adaptativo em comparação com o kernel fixo é grande, ver Tabela 6.5 (G. R. Terrell and Scott 1992b). A minimização foi feita numericamente por Terrell. A vantagem provém da grande porção em forma de sela da densidade.

    Tabela 6.5: Eficiência relativa da estimativa de densidade adaptativa transformada em comparação com o kernel fixo para a densidade normal multivariada.

    A comparação global final é entre o estimador kernel fixo e o estimador de densidade multivariada \(k\)-NN. Este último estimador será tratado como um estimador kernel adaptativo, com o kernel sendo uma densidade uniforme sobre a esfera unitária e a largura de banda sendo considerada como a forma assintótica \[ h= \big( k/(nf(\pmb{x})V_d\big)^{1/d}\cdot \]

    O \(\mbox{AMISE}\) do estimador kernel fixo usando esse mesmo kernel decorre dos resultados gerais de kernel anteriores (ver G. R. Terrell and Scott 1992b). Mostra-se que \[ \dfrac{\mbox{AMISE}_h^*}{\mbox{AAMISE}_k^*}=\left(\dfrac{d-2}{d} \right)^{\frac{d(d+2)}{2(d+4)}} \left(\dfrac{4(d^2-4)}{d^2-6d+16} \right)^{\dfrac{d}{d+4}}\to \dfrac{4}{e}, \] à medida que \(d\to\infty\), ver Tabela 6.6.

    Tabela 6.6: Eficiência relativa do estimador \(k\)-NN em relação ao kernel fixo para a densidade normal multivariada.

    Assim, o estimador de vizinho mais próximo, que se adapta excessivamente às caudas nos casos univariado e bivariado, apresenta melhor desempenho do que o estimador kernel fixo quando \(d\geq 5\), pelo menos para dados normais. Essa superioridade é reconfortante, visto que o algoritmo possui um histórico comprovado em aplicações de alta dimensionalidade, como agrupamento (cluster).


    6.6.4 Algoritmos adaptativos práticos


    Concluímos este tópico com uma discussão sobre dois algoritmos que utilizam a propriedade de não-viesamento do UCV de maneira fundamental. As abordagens PI não conseguem competir com o UCV devido à complexidade assintótica encontrada.


    6.6.4.1 Larguras de banda de biés zero para estimação da cauda


    Uma verdade comum sobre estimativas de densidade não paramétricas é que sua qualidade se degrada nas caudas, onde os dados são escassos. As estimativas geralmente são irregulares em torno de quaisquer pontos de dados que possam ocorrer. No entanto, existe um fenômeno interessante para estimadores de largura de banda fixa que pode fornecer estimativas notavelmente precisas e úteis nas caudas. A apresentação visual também pode ser bastante aprimorada.

    Vamos nos concentrar no comportamento geral do viés de um estimador de kernel positivo quando a largura de banda \(h\) varia de 0 a \(\infty\). Relembre a introdução ao refinamento de dados na Seção 6.6.2.4, onde foi observado que os picos são subestimados e os vales são superestimados.

    Isso é resultado direto do fato de que o viés de uma estimativa kernel é proporcional a \(h^2 f''(x)\). Na Figura 6.27, o viés exato (Fryer 1976) de um estimador kernel normal para uma amostra normal é essencialmente 0 para larguras de banda pequenas, tornando-se negativo ou positivo dependendo se \(0< x< 1\) ou \(x > 1\), respectivamente. Como a densidade \(N(0,1)\) é simétrica, as curvas de viés para valores negativos de \(x\) são idênticas; a densidade é côncava para \(-1 < x < 1\) e convexa para \(|x| > 1\). No outro extremo, quando \(h\to\infty\), a estimativa de densidade \(\widehat{f}(x)\to 0\) para todo \(x\); portanto, o \(\mbox{Viiés}(x)\to -f(x)\), que é negativo em todos os pontos para uma densidade \(N(0,1)\), ver Figura 6.27.

    Figura 6.27: Viés do estimador de densidade kernel normal para uma amostra normal padrão de tamanho \(n = 1.000\) em função da largura de banda.

    Relembrando a equação (6.13), o valor esperado de uma estimativa kernel não depende diretamente do tamanho da amostra \(n\), exceto possivelmente por meio de \(h\). Como a esperança é uma função contínua de \(x\), podemos concluir que, quando \(x > 1\), deve haver uma largura de banda fixa \(h_x\) onde \(\mbox{E}\big(\widehat{f}_{h_x} \big) = f(x)\) exatamente. Essa largura de banda é destacada na Figura 6.27. Referiremos a essa largura de banda como largura de banda de viés zero e a denotaremos por \(h_0(x)\); veja Sain (1994), Sain (2003), Sain and Scott (2002) e também Hazelton (1998), que investigou independentemente o que chamou de largura de banda de aniquilação de viés.

    A existência da largura de banda de viés zero pode resultar em mais de um mínimo local no \(\mbox{MSE}(x)\) pontual e mais de uma largura de banda “ótima” (veja a Figura 6.7). Para \(n\) suficientemente grande, haverá um mínimo local no \(\mbox{MSE}(x)\) próximo a \(h_0(x)\), que denotaremos por \(h^*_0(x)\).

    A largura de banda ótima menor corresponde à teoria assintótica usual, ou seja, \(h^*_x \propto n^{-1/5}\to 0\) quando \(n\to\infty\). No entanto, se optarmos por usar a largura de banda de viés zero \(h_0(x)\), então o viés ao quadrado é 0 e o \(\mbox{MSE}(x)\) é simplesmente a variância \(R(K)f(x)/\big(nh_0(x)\big)\) com boa aproximação, veja as equações (6.11) e (6.15). Essa taxa paramétrica \(O(n^{-1})\) não é alcançável na prática. Tanto Hazelton (1998) quanto Sain (2003) mostram que a taxa é, na verdade, a conhecida \(O(n^{-4/5})\) com dados reais.

    Na Figura 6.28, são mostrados os valores das larguras de banda que (localmente) minimizam o \(\mbox{MSE}(x)\) para a densidade normal padrão e para a densidade de mistura de dois componentes na equação (3.78). No quadro à direita, existe até mesmo uma região próxima a \(x\approx 3.7\) onde existem três dessas larguras de banda.

    Figura 6.28: (Esquerda) Larguras de banda assintoticamente ótimas e larguras de banda de viés zero para uma amostra normal de tamanho \(n = 1000\). Existem duas larguras de banda ótimas quando \(x > 1.218\). A linha tracejada mostra as larguras de banda próximas às de viés zero. (Direita) O mesmo para uma mistura normal, mas com \(n = 500\).

    Na Figura 6.29, comparamos as diversas estimativas kernel disponíveis para uma amostra de tamanho \(n = 500\) da densidade de mistura de dois componentes. O quadro da esquerda utiliza a largura de banda fixa assintoticamente ótima \(h^* = 0.176\), que subsuaviza à esquerda e sobresuaviza à direita.

    Figura 6.29: Diversas estimativas de densidade (linhas contínuas) de uma amostra de tamanho \(n = 500\) a partir da densidade da mistura na equação (3.78), que é mostrada como uma linha tracejada em cada quadro. (Esquerda) Uma estimativa kernel fixa com \(h^* = 0.176\) e a regra de referência normal \(h = 0.473\) (linha pontilhada). (Meio) Estimativa de kernel usando \(h^*(x)\), que é a menor largura de banda no quadro da direita da Figura 6.28. (Direita) Estimativa kernel usando \(h^*_0(x)\) onde existe.

    O quadro do meio na Figura 6.29 mostra o uso de \(h^*(x)\) pontualmente. Isso representa uma melhoria significativa em relação à estimativa kernel fixa, exceto por algumas pequenas anomalias próximas aos pontos de inflexão. Finalmente, no quadro da direita (mesma figura), mostramos as estimativas de viés zero, quando disponíveis. A olho nu, não há erro nessa estimativa, mesmo com uma amostra de apenas 500 pontos. Obviamente, sem o conhecimento das larguras de banda de viés zero exatas, tal resultado não seria esperado. Note que, para o estimador de largura de banda fixa, as larguras de banda da regra de referência normal, UCV, BCV e Sheather-Jones são 0.473, 0.130, 0.222 e 0.218, respectivamente. A largura de banda UCV \(h = 0.130\) (não mostrada) está bem calibrada para a componente direita estreita, deixando a componente esquerda bastante subsuavizada.

    Finalmente, para implementar um procedimento amostral para estimativa de viés zero, a formulação UCV na equação (3.52) pode ser simplesmente modificada incluindo uma função indicadora para dados em um intervalo conhecido assumido \((a,b)\), \[ \int_{-\infty}^\infty \left( \widehat{f}(x)-g(x)\right)^2\pmb{I}_{(a,b)}(x)\mbox{d}x=\int_a^b \widehat{f}(x)^2\mbox{d}x-2\mbox{E}\left(\widehat{f}(X)\pmb{I}_{(a,b)}(X) \right)+c, \] levando a \[ \tag{6.101} \mbox{UCV}_{(a,b)}(h)=\int_a^b \widehat{f}(x)^2 \mbox{d}x-\dfrac{2}{n}\sum_{i=1}^n \widehat{f}_{-i}(x_i)\pmb{I}_{(a,b)}(x_i)\cdot \]

    A implementação mais simples consiste em assumir que a função \(h_0(x)\) pode ser aproximada por uma constante no intervalo \((a,b)\). Obviamente, outras formas funcionais podem ser escolhidas, desde polinômios até splines.

    Retornando ao exemplo da mistura, considere o intervalo \((a,b) = (1.35, 2.51)\) na região convexa em torno do antimodo. O quadro à esquerda na Figura 6.30 exibe a função \(\mbox{UCV}(h)\) amostral na equação (6.101), assumindo uma largura de banda constante em \((a,b)\). Dois mínimos locais são observados, como esperado. A largura de banda maior, \(\widehat{h}^* = 5.71\), está na região de viés zero e poderia ser usada no ponto médio de \((a,b)\) em vez de em todo o intervalo.

    Figura 6.30: (Esquerda) A função \(\mbox{UCV}\) no intervalo \((a,b) = (1.35, 2.51)\) para o conjunto de dados de mistura normal. (Meio) Ajustes log-polinomiais para \(h_0(x)\) que minimizam o \(\mbox{UCV}\) em \((a,b)\). (Direita) A estimativa de viés zero em \((a,b)\) juntamente com a densidade da mistura verdadeira e a estimativa kernel da regra de referência normal.

    Observe, no entanto, que o erro previsto não é menor do que o erro na largura de banda assintoticamente ótima, \(\widehat{h}^* = 0.149\). O quadro do meio exibe os ajustes de largura de banda \(\mbox{UCV}(h)\), onde o logaritmo da largura de banda \(h\) é considerado um polinômio de ordem 0, 1 ou 2 em \((a,b)\).

    Dependendo dos valores iniciais, existem duas soluções, correspondentes aos dois tipos de largura de banda. Para referência, a largura de banda que minimiza o \(\mbox{MSE}(x)\) com viés zero exato, \(h_0(x)\), é mostrada como a linha curva sólida espessa.

    Observe que nenhum dos ajustes polinomiais possui graus de liberdade suficientes para aproximar essa curva. O ajuste de ordem 1 (log-linear) é adequado para \(x < 2.1\), e o ajuste de ordem 2 (log-quadrático) é adequado para \(x > 2.3\). O quadro à direita mostra a estimativa de densidade com viés zero usando o ajuste quadrático, que possui a menor pontuação \(\mbox{UCV}\) 20% menor que uma largura de banda constante. A superestimação da estimativa kernel fixo em \((a,b)\) é corrigida, e a estimativa com viés zero é quase perfeita para \(x > 2\). O ajuste log-linear, não mostrado, também é quase perfeito para \(x < 2\). Um modelo spline para \(h^*_0(x)\) é indicado.

    Concluímos demonstrando uma aproximação simples para a largura de banda de polarização zero, devida a Sain (2003). Como a largura de banda \(h\) é muito grande, utilizamos uma série de Taylor do kernel em torno de 0 na equação (6.12), em vez da mudança de variáveis, e obtemos \[ \begin{array}{rcl} \mbox{E}\big( \widehat{f}(x)\big) & = & \displaystyle \int \dfrac{1}{h}\left(K(0)+\Bigg(\dfrac{x-t}{h}\Bigg)K'(0)+\Bigg(\dfrac{x-t}{h} \Bigg)^2 K''(0)+\cdots \right)f(t)\mbox{d}t \\[0.8em] & = & \displaystyle \dfrac{K(0)}{h}+0+\dfrac{K''(0)}{2h^3}\int \big(x-\mu_f+\mu_f-t \big)^2 f(t)\mbox{d}t+\cdots \\[0.8em] & = & \displaystyle \dfrac{K(0)}{h}+\dfrac{K''(0)}{2h^3}\big((x-\mu_f)^2+\sigma_f^2 \big)+\cdots, \end{array} \] como \(K'(0) = 0\) para kernels simétricos.

    Se \(\mbox{E}\big(\widehat{f}(x)\big) = f(x)\), ou seja, se for não viesado, então \(f(x)\approx K(0)/h_0(x)\) ou \(h_0(x) \approx K(0)/f(x)\), o que é bastante preciso nas caudas da distribuição. Observe que encontrar a função de largura de banda de viés zero está relacionado à estimação da densidade inversa.

    Em resumo, as larguras de banda de viés zero apresentam um potencial interessante para a estimação não paramétrica da densidade em áreas tipicamente consideradas muito difíceis para uma estimação de qualidade. A propriedade de não viés do \(\mbox{UCV}\) permite a investigação dessas larguras de banda excepcionalmente grandes.


    6.6.4.2 Validação Cruzada Unificada (UCV) para estimadores adaptativos


    Considere o estimador de ponto amostral na equação (6.92) e sua óbvia extensão multivariada. Com \(n\) parâmetros, há parâmetros em excesso para validação cruzada prática. Existem diversas estratégias para reduzir o número de parâmetros de suavização, como agrupar os dados em intervalos, nos quais os dados compartilham o mesmo parâmetro de suavização, ou arredondar os dados para os centros dos intervalos. Sain (2002) e Duong and Hazelton (2005) consideram essas opções, bem como parametrizações completas da matriz de suavização \(H\) no caso bivariado.

    Para fins de ilustração, agrupamos os dados não por discretização, mas utilizando o algoritmo \(k\)-means (MacQueen 1967), permitindo que cada cluster tenha seu próprio parâmetro de suavização. O critério usual amostral UCV é otimizado numericamente sobre o vetor de parâmetros de suavização, por exemplo, utilizando a função nlminb do R. A Figura 6.31 exibe as estimativas kernel fixo e adaptativo.

    Figura 6.31: (Esquerda) Doze contornos da estimativa kernel fixo gaussiano bivariado calibrada por UCV \((\widehat{h} = 0.276)\) dos dados padronizados de logaritmo de colesterol e triglicerídeos. (Meio) Sete clusters do \(k\)-means. (Direita) O estimador kernel adaptativo. As sete larguras de banda variam de 0.174 a 2.36. A moda é 54% maior do que no quadro da esquerda. Os 19 níveis de contorno são os mesmos do quadro da esquerda, mais sete em níveis mais altos.

    Sete clusters foram selecionados utilizando \(k\)-means e são mostrados no quadro central. Para reduzir o número de parâmetros, os dados foram padronizados e um kernel gaussiano esférico foi selecionado, o qual possui apenas um parâmetro de suavização (em vez de dois). A largura de banda fixa é 0.174, enquanto as larguras de banda adaptativas são 0.17, 0.25, 0.47, 1.02, 1.48, 1.50 e 2.36. A estimativa adaptativa realça os dois modos principais, bem como o terceiro pico menor, e reduz o ruído nas caudas.

    No entanto, outras escolhas de \(k\) levam a estimativas bastante diferentes, algumas com menos modas, outras com mais. Artefatos também são um problema potencial com estimativas superparametrizadas. Esta foi selecionada por ser mais fiel à estrutura básica do estimador kernel fixo. Ao experimentar com um kernel elíptico totalmente parametrizado, Sain (2002) e Duong and Hazelton (2005) produziram estimativas que poderiam apresentar anomalias ainda mais extremas. Usuários de modelos de mistura normal e do algoritmo de esperança-maximização (EM) estão familiarizados com o potencial de singularidades ao escolher uma matriz de covariância totalmente parametrizada, quando a matriz de covariância ajustada é quase singular. A UCV também é suscetível a tais singularidades.

    Um pacote de mistura bem escrito Mclust, foi aplicado a esses dados (ver Fraley and Raftery 2002). O critério BIC avaliou misturas com \(k\) no intervalo de 1 a 10; \(k = 1\) foi indicado pelo BIC. O ajuste normal bivariado único não detecta a característica bimodal.


    6.7 Aspectos computacionais


    Concluímos nossa análise dos métodos de kernel revisitando diversas opções para calcular um estimador kernel em uma malha igualmente espaçada \(\{a = t_0,t_1,\cdots,t_m = b\}\) com base em uma amostra de tamanho \(n\). Como isso pode ser tornado mais eficiente, como alternativas ao histograma ou ao \(\mbox{ASH}\)?


    6.7.1 Suporte finito do kernel e arredondamento de dados


    A estimativa kernel deve ser avaliada nos pontos médios dos \(m\) intervalos, \(\{m_1, m_2,\cdots, m_m\}\). Se o kernel \(K\), tiver suporte infinito, a abordagem direta requer \(n\cdot m\) avaliações do kernel, ou seja, \(n\) avaliações para cada um dos \(m\) pontos de estimação.

    Entretanto, se um kernel com suporte finito for escolhido, apenas \(f\cdot m\) avaliações do kernel serão necessárias para cada ponto de dados, onde \(f = 2h/(b-a)\), assumindo que o suporte do kernel seja \((-1, 1)\), sem perda de generalidade. Os cálculos devem ser invertidos para que o laço externo percorra os pontos de dados, em vez do laço interno, como mostra este pseudocódigo em R: \[ \delta=(b-a)/m; \quad m_k=\mbox{seq}(a+\delta/2,b-\delta/2,\delta); \quad y=\mbox{rep}(0,m) \\ \mbox{for}(i \mbox{ in } 1:n ) \{k_0=2+[x_i-h-m_1]/\delta; \quad k_1=1+[(x_i+h-m_1)/\delta] \\ \mbox{for}(k \mbox{ in } \max(1,k_0) : \min(m,k_1))\{ y_k=y_k+K_h(x_i-m_k)\}\}; \quad y=y/n\cdot \]

    No entanto, o trabalho ainda é proporcional ao tamanho da amostra \(n\). Pré-agrupar ou arredondar os dados para a mesma malha reduz o trabalho de \(f\cdot m\cdot n\) para \(f\cdot m^2\).

    A parte do código que precisa ser modificada é: \[ \mbox{for}(\ell \mbox{ in } 1:m ) \{k_0=\ell+1-[h/\delta]; \quad k_1=\ell+[h/\delta] \\ \mbox{for}(k \mbox{ in } \max(1,k_0) : \min(m,k_1))\{ y_k=y_k+K_h(m_\ell-m_k)\}\}\cdot \]

    O uso de dados arredondados em um estimador kernel de densidade atraiu muita atenção em relação à perda de precisão, ver Scott (1981), Scott and Sheather (1985) e Silverman (1982); ver também Jones (1989). Härdle and Scott (1988) cunharam a expressão “média ponderada de pontos arredondados” ou método WARPing para descrever a abordagem geral e aplicaram o método WARPing a outros algoritmos multivariados.

    Com os computadores modernos, \(m\) pode ser bastante grande e \(\delta\) bastante pequeno, de modo que o erro de aproximação introduzido seja desprezível. Em múltiplas dimensões, no entanto, manter uma grade inteira no núcleo deve ser considerado quando \(d\geq 3\).


    6.7.2 Convolução e transformadas de Fourier


    Na Seção 6.1.2, exploramos uma motivação para o estimador de kernel fixo via convolução. Essa abordagem de engenharia também leva a algoritmos computacionais eficientes que utilizam métodos de Fourier para realizar a convolução.

    Começamos com três definições: \[ \tag{6.102} \mbox{Convolução: } \qquad (f*g)(x)=\int_{-\infty}^\infty f(t)g(x-t)\mbox{d}t, \] \[ \tag{6.103} \mbox{Fourier Xfm: } \qquad \widetilde{f}(\xi) = \int_{-\infty}^\infty f(x)\exp\big(-2\pi \, i \, x \,\xi \big)\mbox{d}x, \] e \[ \tag{6.104} \mbox{Fourier Xfm Inversa: } \qquad f(x)=\int_{-\infty}^\infty \widetilde{f}(\xi)\exp\big(2\pi \, i \, \xi \,x \big)\mbox{d}\xi\cdot \]

    Por exemplo, representando os dados \(\{x_1,\cdots,x_n\}\) através da função de ddensidade empírica \(f_n(x)\), em (2.2), \[ \tag{6.105} \begin{array}{rcl} \widetilde{f}_n(\xi) & = & \displaystyle \int \left(\dfrac{1}{n}\sum_{j=1}^n \delta(x-x_j) \right) \exp\big(-2\pi \, i \, x \,\xi \big)\mbox{d}x\\[0.8em] & =& \displaystyle \dfrac{1}{n}\sum_{j=1}^n \int \delta(x-x_j)\exp\big(-2\pi \, i \, x \,\xi \big) \mbox{d}x=\dfrac{1}{n}\sum_{j=1}^n \exp\big(-2\pi \, i \, x_j \,\xi \big)\cdot \end{array} \]

    Uma alternativa à realização da operação de convolução direta (6.102), que pode ser muito lenta, é realizar três Transformadas de Fourier para obter \(c(x)=(f*g)(x)\) por: \[ \widetilde{c}(\xi)=\widetilde{f}(\xi)\cdot \widetilde{g}(\xi) \] o qual implica em \[ c(x)=\widetilde{\widetilde{c}(\xi)}=\widetilde{\widetilde{f}(\xi)\cdot \widetilde{g}(\xi)}\cdot \]


    6.7.2.1 Aplicação a estimadores de densidade kernel


    Aplicada a um estimador kernel, a equação (6.102) torna-se \[ \tag{6.106} \widetilde{\widehat{f}}(x)=\int \widehat{f}(x)\exp\big(-2\pi \, i \, x \,\xi \big)\mbox{d}x, \] a qual calcula-se como \[ \begin{array}{rcl} \widetilde{\widehat{f}}(x) & = & \displaystyle \int \left(\dfrac{1}{n}\sum_{j=1}^n \dfrac{1}{h}K\Bigg(\dfrac{x-x_j}{h} \Bigg) \right) \exp\big(-2\pi \, i \, x \,\xi \big)\mbox{d}x \\[0.8em] & = & \displaystyle \dfrac{1}{n}\sum_{j=1}^n \left(\int K(s_j) \exp\big(-2\pi \, i \,(x_j+hs_j) \,\xi \big) \mbox{d}s_j\right) \\[0.8em] & = & \displaystyle \dfrac{1}{n}\sum_{j=1}^n \left( \exp\big(-2\pi \, i \, x_j \,\xi \big) \int K(s_j) \exp\big(-2\pi \, i \,(hs_j) \,\xi \big) \mbox{d}s_j\right) \\[0.8em] \end{array} \] que escreve-se como \[ \tag{6.107} \widetilde{\widehat{f}}(x) = \dfrac{1}{n}\sum_{j=1}^n \exp\big(-2\pi \, i \, x_j \, \xi \big) \widetilde{K}(hs_j) \] ou \[ \widetilde{\widehat{f}}(x) = \widetilde{f}_n \cdot \widetilde{K}(hs), \] utilizando a mudança de variáveis \(s_j = (x-x_j)/h\) e observando que \(\widetilde{K}(hs_j)\) é o mesmo para todo \(j\).

    Agora, \(\widetilde{K}\) é conhecido para muitos kernels, incluindo o kernel gaussiano: \[ \widetilde{\phi}_h(\xi) = \exp\big(-2\pi^2 \, h^2 \, \xi^2 \big)\cdot \] Portanto, (6.107) torna-se \[ \tag{6.108} \widetilde{\widehat{f}}(\xi) = \dfrac{1}{n}\sum_{j=1}^n \exp\big(-2\pi \, i \, x_j \, \xi \big) \exp\big(-2\pi^2 \, h^2 \, \xi^2 \big)\cdot \]

    Podemos verificar a correção de (6.108) calculando sua Transformação Inversa de Fourier usando o Mathematica (2012): \[ \tag{6.109} \widetilde{\widetilde{\widehat{f}}}(x) = \widehat{f}(x) = \int \widetilde{\widehat{f}}(\xi)\exp\big(2\pi \, i \, x \, \xi \big) \mbox{d}\xi \] a qual escreve-se como \[ \tag{6.110} \widetilde{\widetilde{\widehat{f}}}(x) = \dfrac{1}{n}\sum_{j=1}^n \int \exp\big(-2\pi \, i \, x_j \, \xi \big) \exp\big(-2\pi^2 \, h^2 \, \xi^2 \big)\exp\big(2\pi \, i \, x \, \xi \big)\mbox{d}\xi \] ou \[ \tag{6.111} \widetilde{\widetilde{\widehat{f}}}(x) = \dfrac{1}{n}\sum_{j=1}^n \dfrac{1}{\sqrt{2\pi}h} \exp\big(-(x-x_j)^2/2h^2\big) \] que é exatamente a estimativa da densidade kernel gaussiano.


    6.7.2.2 FFTs


    Na prática, a Transformada Discreta de Fourier e a FFT são empregadas para calcular uma aproximação muito boa de \(\widehat{f}(x)\).

    A FFT é aplicada a um intervalo finito, digamos \(\big(-\frac{1}{2},\frac{1}{2}\big)\) sem perda de generalidade e a estimativa resultante é periódica. Para evitar efeitos de borda, os dados são reescalados para um subintervalo, como \((-0.3, 0.3)\). Silverman (1986) recomenda um buffer de pelo menos \(3h\) em cada extremidade.

    No contexto de séries temporais, o preenchimento com zeros na estimativa da densidade espectral é um requisito bem conhecido para evitar aliasing, ou seja, sobreposição de frequências (Blackman and Tukey 1958).

    O agrupamento dos dados é fundamental para a abordagem, ou seja, \[ \widehat{f}(x)=\dfrac{1}{n}\sum_{i=1}^n K_h(x-x_i)\approx \dfrac{1}{n}\sum_{k=1}^m \nu_k K_h(x-m_k), \] e se calcularmos \(\widehat{f}(x)\) nos mesmos pontos médios dos intervalos \[ \tag{6.112} \widehat{f}(m_\ell)=\dfrac{1}{n}\sum_{k=1}^m \nu_k K_h(m_\ell-m_k)=\sum_{k=1}^m \nu_k \dfrac{1}{n}K_h\big(\delta(\ell-k)\big), \] que é precisamente uma convolução discreta.

    Em R, isso pode ser implementado através da função convolve, com argumentos fornecidos pelas contagens de intervalos e os valores do kernel em \(t = m_\ell\) divididos por \(n\). Para que os valores finais estejam devidamente alinhados, os valores do kernel devem ser deslocados circularmente pela metade.

    A extensão multivariada da abordagem FFT é discutida por Wand (1994) e Wand and Jones (1995). Esses autores e Silverman (1986) defendem a substituição do agrupamneto simples pelo agrupamento linear, embora, com a computação moderna e malhas muito finas, o impacto possa ser bastante pequeno.


    6.7.2.3 Discussão


    A função de densidade padrão no R agora utiliza o algoritmo FFT no caso univariado. Em situações onde um kernel de fronteira é necessário, outros algoritmos precisarão ser empregados. O uso da FFT em mais de uma dimensão também é possível, mas a malha cresce exponencialmente, de modo que três dimensões representam um limite prático.

    A abordagem \(\mbox{ASH}\) vai além de três dimensões, calculando fatias ou densidades condicionais. Também é possível analisar subconjuntos do domínio. De qualquer forma, existem muitas outras oportunidades para experimentar kernels de estimadores \(\mbox{ASH}\) no contexto multivariado, e o leitor é encorajado a testá-los em seus dados. Além dos softwares \(\mbox{ASH}\) e ashn disponíveis no R e fornecidos pelo autor, o pacote ks é recomendado (Duong 2007). Este pacote utiliza a biblioteca de animação rgl.

    Este capítulo abordou os métodos kernel, que muitas vezes constituem o único tópico de um livro. A abordagem kernel foi motivada de diferentes maneiras que podem interessar a um estatístico, analista numérico, engenheiro ou matemático. Questões relacionadas à escolha do kernel — sua ordem, propriedades de fronteira e projeto — foram analisadas. Extensões para dados multivariados e o uso do kernel produto foram revisados.

    Ideias de validação cruzada discutidas em capítulos anteriores foram estendidas aos kernels. Finalmente, o importante tópico de métodos adaptativos locais foi introduzido. Os métodos adaptativos são muito promissores, mas geralmente introduzem muitos parâmetros novos que são difíceis de estimar e, frequentemente, introduzem artefatos da amostra, em vez da densidade subjacente.

    Uma interessante revisão é fornecida por Jones and Signorini (1997), que abordam esses algoritmos adaptativos e outros no contexto de estratégias de viés de ordem superior. Eles observam que há muitos detalhes técnicos a serem superados. Examinamos apenas alguns. Uma estratégia conservadora que pode ser recomendada é manter uma estimativa kernel de largura de banda fixa e estar ciente e atento às suas limitações.


    6.8 Exercícios


    • 1- Calcule o \(\mbox{IV}\) do estimador (6.2) e compare com o resultado do histograma. Qual é o kernel se \[ \widehat{f}(x)=\dfrac{1}{h}\big(F_n(x+h)-F_n(x) \big) \] for utilizado?


    • 2- Calcule o viés e a variância diretamente para o estimador (6.3). Demonstre isso com o kernel equivalente \(K=U(-0.5,0.5)\).


    • 3- Verifique se os resultados \(\mbox{ASH}\) são obtidos a partir do Teorema 6.1 quando o kernel do triângulo isósceles é especificado.


    • 4- Examine os gráficos do kernel equivalente de Wahba (6.9) para várias escolhas dos dois parâmetros \(p\) e \(\lambda\).


    • 5- Após a discussão do rootgrama para o histograma, mostre que a raiz quadrada da estimativa kernel estabiliza a variância em comparação com (6.15). Mostre que a variância da estimativa kernel da raiz é \(R(K)/\big(4nh\big)\).


    • 6- Verifique as equações (6.21). Encontre as larguras de banda ótimas para estimar a primeira e a segunda derivadas de \(f\). Avalie-as para a densidade \(N(\mu,\sigma^2)\).


    • 7- Considere o kernel \(K_4(t)\) na Tabela 6.1. E se você suavizar esse kernel aumentando a potência do primeiro fator, como \[ c(1-t^2)^k (a+bt^2), \] onde \(k\) é 2 ou 3, e as constantes \(c\), \(a\) e \(b\) são escolhidas de forma que o kernel seja de ordem 4 e integre 1? Como esses kernels alternativos se comportam teoricamente? E visualmente?


    • 8- Mostre que duas estimativas kernel com tamanhos de amostra na proporção dada na equação (6.26) têm o mesmo \(\mbox{AMISE}\).


    • 9- Mostre que o \(\mbox{IV}\) da estimativa kernel fixo é igual (exatamente) \[ \mbox{IV}(h)=\dfrac{1}{h}\int \mbox{E}\big(K_h(x-X)^2 \big)\mbox{d}x - \dfrac{1}{n}\int \mbox{E}^2\big(K_h(x-X) \big)\mbox{d}x\cdot \] Mostre que o primeiro termo no \(\mbox{IV}\) é exatamente \(R(K)/(nh)\). Mostre que os termos seguintes são \[ -\dfrac{R(f)}{n}+\dfrac{h^2 \mu_2 R(f')}{n}-\dfrac{h^4}{n}\Bigg(\dfrac{\mu_2^2}{4}+\dfrac{\mu_4}{24} \Bigg)R(f'')+\cdots, \] onde \(\mu_k\) é o \(k\)-ésimo momento do kernel.


    • 10- Mostre que o \(\mbox{ISB}\) de uma estimativa kernel fixa é \[ \mbox{ISB}(h) = \int \Big( \mbox{E}\big(K_h(x-X) \big) -f(x)\Big)^2\mbox{d}x \cdot \] Suponha que o kernel \(K\) seja simétrico em torno de zero, de modo que \(\displaystyle \int \omega^k K(\omega)\mbox{d}\omega = 0\) para \(k\) ímpar. Mostre que os primeiros termos de viés são iguais a \[ h^4 \dfrac{\mu_2^2}{4}R(f'')-h^6 \dfrac{\mu_2 \mu_4}{24}R(f''')+h^8 \Bigg(\dfrac{\mu_4^2}{476}+\dfrac{\mu_2\mu_6}{720} \Bigg) R(f^{vi})\cdot \] Dica: Observe os termos do produto cruzado e note que \(\displaystyle \int f'' f^{iv}=-\int (f''')^2\), por exemplo.


    • 11- Verifique as expressões de viés para os estimadores de diferenças finitas de ordem superior na equação (6.32). Elabore seus próprios kernels boxcar de ordem superior usando espaçamentos diferentes de múltiplos inteiros de \(h\).


    • 12- Compare os valores de \(\mbox{AMSE}(x)\) para quatro combinações de estimativas pontuais do kernel, em 0 e 1 com um kernel de segunda e quarta ordem.


    • 13- Compare empiricamente o método kernel de quarta ordem com o estimador de razão de quarta ordem de Terrell-Scott para dados normais simulados, bem como para os dados de queda de neve (snow).


    • 14- Derive o estimador de razão kernel de quarta ordem de Terrell-Scott. Estenda o procedimento para um estimador de sexta ordem.


    • 15- Calcule o kernel equivalente do estimador paramétrico \(\phi(x \,| \,0, s^2)\).


    • 16- Calcule o \(k\)-ésimo momento “teórico” de \(\widehat{f}\), tratando o estimador kernel como uma densidade “verdadeira”.


    • 17- Encontre o kernel do polígono de frequência indiferente mencionado na Tabela 6.2.


    • 18- Considere a classe de kernels Beta deslocados \(c_k (1-t^2)^k_+\). Encontre sua variância e reescale-os de forma que cada um tenha variância 1. Mostre que esses kernels reescalados convergem para um kernel normal padrão quando \(k\to\infty\).


    • 19- Derive os resultados \(\mbox{AMISE}\) do kernel produto a partir das fórmulas gerais do kernel multivariado.


    • 20- Verifique a fórmula da largura de banda equivalente para kernels de ordem superior na equação (6.31). Verifique se o último fator é aproximadamente igual a 1 para vários kernels.


    • 21- Usando uma série de Taylor em \(\Delta F_n(x, kh)\), verifique se as estimativas de diferenças finitas na equação (6.32) são de ordem superior. Derive algumas estimativas adicionais com base nos espaçamentos \(h\), \(2h\), \(4h\), \(8h\),\(\cdots\)


    • 22- Verifique por integração direta se a “variância” do kernel na equação (6.36) é zero.


    • 23- Encontre kernels de modificação de fronteira com base no kernel de Epanechnikov. Compare com kernels suportados em \([c, 1]\). Investigue o aumento em \(R(K_c)\). Qual deve ser a largura (à direita) do kernel para que a rugosidade seja a mesma do kernel biweight? Existe sempre uma “rugosidade equivalente” desse tipo?


    • 24- Verifique a equação (6.40).


    • 25- Relembrando o estimador (6.45). Mostre que este é funcionalmente equivalente a escolher \(K\) como \(N(\pmb{0}_d,\Sigma)\) com \(H = \pmb{I}_d\), ou escolher \(K\) como \(N(\pmb{0}_d,\pmb{I}_d)\) com \(H = \Sigma^{1/2}\). Assim, a transformação linear pode ser aplicada tanto aos dados quanto ao kernel, conforme a preferência.


    • 26- Kernels com suporte finito não precisam ter apenas um número finito de derivadas. Por exemplo, considere \[ K(t)\propto e^{-1/(1-t^2)}\pmb{I}_{(-1,1)}(t)\cdot \] Mostre que a constante de normalização é \(\sqrt{\pi}/e\times \mbox{Hipergeométrica U} \Big[\frac{1}{2},0,1 \Big]\). Trace o gráfico do kernel. Dica: Use a mudança de variável \(t = (1-x^2)^{-1}-1\) para \(x\in (0, 1)\).


    • 27- Mostre que \(c \, (x-a)^4/144\) é a solução da equação diferencial da densidade adaptativa nula (6.90).


    • 28- Verifique o kernel equivalente para o estimador paramétrico \(N(\overline{x}, 1)\). Trace o gráfico de \(K(x,t)\). Calcule o kernel equivalente para o estimador paramétrico \(N(0, s^2 )\) e trace o gráfico.


    • 29- Lembre-se de que \(\sigma^2_{\widehat{f}} = s^2_x + h^2 \sigma_K^2\). A partir desse resultado, argumente que se \(s_x > \sigma_f\), então \(\widehat{h}_{ISE} > \widehat{h}_{MISE}\) é o resultado provável. O inverso é verdadeiro?


    • 30- Experimente o estimador de séries ortogonais simples em alguns dados Beta com \(0\leq m\leq 6\). Com uma amostra maior, há controle suficiente na estimativa apenas com \(m\)?


    • 31- Usando a estimativa simples para os coeficientes de Fourier na equação (6.6), encontre estimativas não viesadas dos dois termos desconhecidos na equação (6.62).


    • 32- Mostre que se um fator de correção da forma \((1+b\, n^{-\delta})\) for aplicado a \(h^*\) em (6.79), então a melhor escolha é \(\delta = 1/5\).


    • 33- Deduza a equação (6.85). Dica: A fração de pontos, \(k/n\), na bola de raio \(h\) centrada em \(\pmb{x}\) é aproximadamente igual a \(f(\pmb{x})\) vezes o volume da bola.


    • 34- Mostre que o estimador \(k\)-NN usando a função de distância \(d_k\) tem integral infinita. Dica: No caso univariado, observe o estimador para \(x> x_{(n)}\), a maior estatística de ordem.


    • 35- Use a desigualdade de Jensen para mostrar que a razão na equação (6.89) é \(\leq 1\). Dica: A integral entre parênteses no denominador é \[ \mbox{E}\big(f^{(p)}(X)^2\big/f(X) \big)\cdot \]


    • 36- (Pesquisa). Em vez de abandonar o \(\mbox{UCV}\) em favor de estimadores assintóticos mais eficientes, considere ajustar os dados antes de calcular a curva \(\mbox{UCV}\). O ruído no \(\mbox{UCV}\) deve-se, em parte, ao fato de os termos \((x_i-x_j)/h\) serem muito ruidosos para \(|x_i-x_j| < h\). Imagine colocar pequenas molas entre os pontos de dados, resultando em uma variação muito suave das distâncias entre os pontos. O \(\mbox{UCV}\) resultante parece ter um comportamento muito melhor. Um método adaptativo local substitui \(x_i\gets (x_{i-1}+x_{i+1})/2\), assumindo que os dados foram ordenados. Experimente essa modificação em dados simulados.


    • 37- Prove que o viés médio de uma estimativa kernel é exatamente 0. Dica: O viés pontual é \[ \int K(\omega) f(x-h\omega)\mbox{d}\omega - f(x)\cdot \]


    • 38- Wand and Jones (1995) usam duas identidades matriciais, que podem ser encontradas em Neudecker and Magnus (1988), para derivar a equação (6.51). Suponha que \(B\) e \(C\) sejam matrizes simétricas \(d\times d\). Então \(\mbox{hvec}(B)\) vetoriza a porção triangular inferior, enquanto \(\mbox{vec}(B)\) vetoriza a matriz inteira. Esses dois vetores são conectados pela matriz de duplicação \(D_d\) de dimensão \(d^2\times \frac{d(d+1)}{2}\), cujas entradas são todas 0 e 1; ou seja, \[ D_d \, \mbox{vech}(B) = \mbox{vec}(B)\cdot \] As identidades e relações são \[ \begin{array}{rcl} \mbox{vec}(B) & = & D_d \, \mbox{vech}(B),\\[0.8em] D^\top_d \, \mbox{vec}(B) & = & \mbox{vech}\big( 2B-\mbox{Diag}(B)\big),\\ \mbox{tr}(B^\top C) & = & \mbox{vec}(B)^\top\mbox{vec}(C) = \mbox{vec}(C)^\top \mbox{vec}(B)\cdot \end{array} \] Complete o seguinte argumento. Seja \(B = AA^\top\) na equação (6.50). Escreva o integrando como \[ \begin{array}{rcl} \mbox{tr}\big(B \, \nabla^2 f(\pmb{x})\big)^2 & = & \displaystyle \mbox{vec}(B)^\top \big(\mbox{vec}(\nabla^2 f(\pmb{x})) \big)\times \big(\mbox{vec}(\nabla^2 f(\pmb{x})) \big)^\top \mbox{vec}(B) \\[0.8em] & = & \displaystyle \mbox{vech}(B)^\top D^\top_d \big(\mbox{vec}(\nabla^2 f(\pmb{x})) \big)\times \big(\mbox{vec}(\nabla^2 f(\pmb{x})) \big)^\top D_d \, \mbox{vech}(B)\cdot \end{array} \] A dedução pode ser concluída usando a segunda identidade em \(D^\top_d\big(\mbox{vec}(\nabla ^2 f(\pmb{x})) \big)\).

    6.9 Bibliografia


    Abramson, I. S. 1982. “On Bandwidth Variation in Kernel Estimates—a Square Root Law.” Annals of Statistics, no. 10: 1217–23.
    Altman, N. S. 1990. “Kernel Smoothing of Data with Correlated Errors.” Journal of the American Statistical Association, no. 85: 749–59.
    Bartlett, M. S. 1963. “Statistical Estimation of Density Functions.” Sankhyā, Ser. A, no. 25: 245–54.
    Blackman, R. B., and J. W. Tukey. 1958. Measurement of Power Spectra.
    Bowman, A. W. 1984. “An Alternative Method of Cross-Validation for the Smoothing of Density Estimates.” Biometrika, no. 71: 353–60.
    ———. 1985. “A Comparative Study of Some Kernel-Based Nonparametric Density Estimators.” Journal of Statistical Computation and Simulation, no. 21: 313–27.
    Bowyer, A. 1980. “Experiments and Computer Modelling in Stick-Slip.” PhD thesis, University of London, England.
    Box, G. E. P., and D. R. Cox. 1964. “An Analysis of Transformations.” Journal of the Royal Statistical Society, B, no. 26: 211–43.
    Breiman, L., W. Meisel, and E. Purcell. 1977. “Variable Kernel Estimates of Multivariate Densities and Their Calibration.” Technometrics, no. 19: 135–44.
    Cacoullos, T. 1966. “Estimation of a Multivariate Density.” Annals of the Institute of Statistical Mathematics, no. 18: 178–89.
    Cencov, N. N. 1962. “Evaluation of an Unknown Density from Observations.” Soviet Mathematics, no. 3: 1559–62.
    Chiu, S. T. 1989. “Bandwidth Selection for Kernel Estimates with Correlated Noise.” Statistical & Probility Letters, no. 8: 347–54.
    Choi, E., and P. Hall. 1999. “Data Sharpening as a Prelude to Density Estimation.” Biometrika, no. 86: 941–47.
    Chow, Y. S., S. Geman, and L. D. Wu. 1983. “Consistent Cross-Validated Density Estimation.” Annals of Statistics, no. 11: 25–38.
    Davis, K. B. 1975. “Mean Square Error Properties of Density Estimates.” Annals of Statistics, no. 3: 1025–30.
    Deheuvels, P. 1977a. “Estimation Non Paramétrique de La Densité Par Histogrammes Généralises.” Revue de Statistique Appliquée, v.15, no. 25/3: 5–42.
    ———. 1977b. “Estimation Non Paramétrique de La Densité Par Histogrammes Généralises II.” Publications de l’Institute Statistique de l’Université Paris, no. 22: 1–23.
    Duin, R. P. W. 1976. “On the Choice of Smoothing Parameter for Parzen Estimators of Probability Density Functions.” IEEE Transactions on Computers, no. C-25: 1175–79.
    Duong, T. 2007. “Ks: Kernel Density Estimation and Kernel Discriminant Analysis for Multivariate Data in r.” Journal of Statistical Software, no. 21(7): 1–16.
    Duong, T., and M. L. Hazelton. 2005. “Convergence Rates for Unconstrained Bandwidth Matrix Selectors in Multivariate Kernel Density Estimation.” Journal of the Multivariate Analysis, no. 93: 417–33.
    Epanechnikov, V. K. 1969. “Non-Parametric Estimation of a Multivariate Probability Density.” Theory of Probability and Its Applications, no. 14: 153–58.
    Farrell, R. H. 1972. “On the Best Obtainable Asymptotic Rates of Convergence in Estimation of the Density Function at a Point.” Annals of Mathematical Statistics, no. 43: 170–80.
    Fix, E., and J. L.Jr. Hodges. 1951. “Nonparametric Discrimination: Consistency Properties.” USAF School of Aviation Medicine, Randolph Field, Texas.
    Fraley, C., and A. E. Raftery. 2002. “Model-Based Clustering, Discriminant Analysis, and Density Estimation.” Journal of the American Statistical Association, no. 97: 611–31.
    Fryer, M. J. 1976. “Some Errors Associated with the Non-Parametric Estimation of Density Functions.” Journal of the Institute of Mathematics and Its Applications, no. 18: 371–80.
    Gasser, T., and H. G. Müller. 1979. “Kernel Estimation of Regression Functions.” Springer-Verlag, Berlin.
    Gasser, T., H. G. Müller, and V. Mammitzsch. 1985. “Kernels for Nonparametric Curve Estimation.” Journal of the Royal Statistical Society, B, no. 47: 238–52.
    Good, I. J., and R. A. Gaskins. 1972. “Global Nonparametric Estimation of Probability Densities.” Virginia Journal of Science, no. 23: 171–93.
    Habbema, J. D. F., J. Hermans, and K. Van Der Broek. 1974. “A Stepwise Discriminant Analysis Program Using Density Estimation.” Physica-Verlag, Vienna.
    Hall, P. 1984. “Central Limit Theorem for Integrated Square Error of Multivariate Density Estimators.” Journal of the Multivariate Analysis, no. 14: 1–16.
    Hall, P., and J. S. Marron. 1987a. “Extent to Which Least-Squares Cross-Validation Minimises Integrated Square Error in Nonparametric Density Estimation.” Probability Theory and Related Fields, no. 74: 567–81.
    ———. 1987b. “On the Amount of Noise Inherent in Bandwidth Selection for a Kernel Density Estimator.” Annals of Statistics, no. 15: 163–81.
    ———. 1988. “Variable Window Width Kernel Estimates of Probability Densities.” Probability Theory and Related Fields, no. 80: 37–49.
    ———. 1991. “Lower Bounds for Bandwidth Selection in Density Estimation.” Probability Theory and Related Fields, no. 90: 149–73.
    Hall, P., and M. C. Minnotte. 2002. “High Order Data Sharpening for Density Estimation.” Journal of the Royal Statistical Society, B.
    Härdle, W., and D. W. Scott. 1988. “Smoothing in Low and High Dimensions by Weighted Averaging Using Rounded Points.” Computational Statistics, no. 7: 97–128.
    Hazelton, M. L. 1998. “Bias Annihilating Bandwidths for Kernel Density Estimation at a Point.” Statistics & Probability Letters, no. 38: 305–9.
    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.
    Jones, M. C. 1989. “Discretized and Interpolated Kernel Density Estimates.” Journal of the American Statistical Association, no. 84: 733–41.
    ———. 1990. “Variable Kernel Density Estimates.” Australian Journal of Statistics, no. 32: 361–71.
    Jones, M. C., and R. F. Kappenman. 1992. “On a Class of Kernel Density Estimate Bandwidth Selectors.” Scandinavian Journal of Statistics, no. 19: 337–49.
    Jones, M. C., and D. F. Signorini. 1997. “A Comparison of Higher-Order Bias Kernel Density Estimators.” Journal of the American Statistical Association, no. 92: 1064–73.
    Klonias, V. K. 1982. “Consistency of Two Nonparametric Maximum Penalized Likelihood Estimators of the Probability Density Function.” Annals of Statistics, no. 10: 811–24.
    Kronmal, R. A., and M. E. Tarter. 1968. “The Estimation of Probability Densities and Cumulatives by Fourier Series Methods.” Journal of the American Statistical Association, no. 63: 925–52.
    Loftsgaarden, D. O., and C. P. Quesenberry. 1965. “A Nonparametric Estimate of a Multivariate Density Function.” Annals of the Mathematical Statistics, no. 36: 1049–51.
    MacQueen, J. B. 1967. “Some Methods for Classification and Analysis of Multivariate Observations.” University of California Press.
    Marron, J. S., and D. Nolan. 1988. “Canonical Kernels for Density Estimation.” Statistics and Probability Letters, no. 7: 195–99.
    Marron, J. S., and M. P. Wand. 1992. “Exact Mean Integrated Squared Error.” Annals of Statistics, no. 20: 712–36.
    Montricher, G. F., R. A. Tapia, and J. R. Thompson. 1975. “Nonparametric Maximum Likelihood Estimation of Probability Densities by Penalty Function Methods.” Annals of Statistics, no. 3: 1329–48.
    Müller, H. G. 1988. Nonparametric Regression Analysis of Longitudinal Data. Springer-Verlag, Berlin.
    Neudecker, H., and J. R. Magnus. 1988. Matrix Differential Calculus with Applications in Statistics and Econometrics. John Wiley, New York.
    Nezames, D. 1980. “Some Results for Estimating Bivariate Densities Using Kernel, Orthogonal Series, and Penalized-Likelihood Procedures.” PhD thesis, Department of Mathematical Sciences, Rice University.
    Park, B. U., and J. S. Marron. 1990. “Comparison of Data-Driven Bandwidth Selectors.” Journal of the America Statistical Association, no. 85: 66–72.
    Parzen, E. 1962. “On Estimation of Probability Density Function and Mode.” Annals of the Mathematical Statistics, no. 33: 1065–76.
    Rice, J. A. 1984. “Boundary Modification for Kernel Regression.” Communications in Statistics, no. 13: 893–900.
    Rosenblatt, M. 1956. “Remarks on Some Nonparametric Estimates of a Density Function.” Annals of Mathematical Statistics, no. 27: 832–37.
    Rudemo, M. 1982. “Empirical Choice of Histograms and Kernel Density Estimators.” Scandinavian Journal of Statistics, no. 9: 65–78.
    Sain, S. R. 1994. “Adaptive Kernel Density Estimation.” PhD thesis, Department of Statistics; Rice University.
    ———. 2002. “Multivariate Locally Adaptive Density Estimation.” Computational Statistics & Data Analysis, no. 39: 165–86.
    ———. 2003. “A New Characterization and Estimation of the Zero–Bias Bandwidth.” The Australian & New Zealand Journal of Statistics, no. 45: 29–42.
    Sain, S. R., K. A. Baggerly, and D. W. Scott. 1994. “Cross-Validation of Multivariate Densities.” Journal of the American Statistical Association, no. 89: 807–17.
    Sain, S. R., and D. W. Scott. 2002. “Zero-Bias Locally Adaptive Density Estimators.” Scandinavian Journal of Statistics, no. 29: 441–60.
    Samiuddin, M., and G. M. El-Sayyad. 1990. “On Nonparametric Kernel Density Estimates.” Biometrika, no. 77: 865–74.
    Schucany, W. R. 1989. “Locally Optimal Window Widths for Kernel Density Estimation with Large Samples.” Statistics and Probability Letters, no. 7: 401–5.
    Schucany, W. R., and J. P. Sommers. 1977. “Improvement of Kernel Density Estimators.” Journal of the American Statistical Association, no. 72: 420–23.
    Schuster, E. F., and G. G. Gregory. 1981. “On the Nonconsistency of Maximum Likelihood Nonparametric Density Estimators.” Springer-Verlag, New York.
    Schwartz, S. C. 1967. “Estimation of a Probability Density by an Orthogonal Series.” Annals of the Mathematical Statistics, no. 38: 1262–65.
    Scott, D. W. 1976. “Nonparametric Probability Density Estimation by Optimization Theoretic Techniques.” PhD thesis, Department of Mathematical Sciences, Rice University.
    ———. 1981. “Using Computer-Binned Data for Density Estimation.” Springer-Verlag, New York.
    ———. 1985. “Averaged Shifted Histograms: Effective Nonparametric Density Estimators in Several Dimensions.” Annals of Statistics, no. 13: 1024–40.
    ———. 1988. “Comment on w. Härdle, p. Hall, and j.s. Marron “How Far Are Automatically Chosen Regression Smoothing Parameters from Their Optimum.” Journal of the American Statistical Association, no. 83: 96–98.
    ———. 2015. Multivariate Density Estimation. John Wiley & Sons. Inc.
    Scott, D. W., and L. E. Factor. 1981. “Monte Carlo Study of the Three Data-Based Nonparametric Density Estimators.” Journal of the American Statistical Association, no. 76: 9–15.
    Scott, D. W., A. M. Gotto, J. S. Cole, and G. A. Gorry. 1978. “Plasma Lipids as Collateral Risk Factors in Coronary Artery Disease: A Study of 371 Males with Chest Pain.” Journal of Chronic Diseases, no. 31: 337–45.
    Scott, D. W., and S. J. Sheather. 1985. “Kernel Density Estimation with Binned Data.” Communications in Statistics, no. 14: 1353–59.
    Scott, D. W., R. A. Tapia, and J. R. Thompson. 1977. “Kernel Density Estimation Revisited.” Nonlinear Analysis: Theory, Methods & Applications, no. 1: 339–72.
    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.
    Sheather, S. J. 1983. “A Data-Based Algorithm for Choosing the Window Width When Estimating the Density at a Point.” Computational Statistics & Data Analysis, no. 1: 229–38.
    Sheather, S. J., and M. C. Jones. 1991. “A Reliable Data-Based Bandwidth Selection Method for Kernel Density Estimation.” Journal of the Royal Statistical Society, B, no. 53: 683–90.
    Silverman, B. W. 1978a. “Choosing the Window Width When Estimating a Density.” Biometrika, no. 65: 1–11.
    ———. 1978b. “Density Ratios, Empirical Likelihood and Cot Death.” Journal of Applied Statistics, no. 27: 26–33.
    ———. 1982. “Algorithm AS176. Kernel Density Estimation Using the Fast Fourier Transform.” Journal of Applied Statistics, no. 31: 93–99.
    ———. 1986. Density Estimation for Statistics and Data Analysis. Chapman; Hall, London.
    Silverman, B. W., and M. C. Jones. 1989. “E. Fix and JL Hodges (1951): An Important Contribution to Nonparametric Discriminant Analysis and Density Estimation: Commentary on Fix and Hodges (1951).” International Statistical Review, no. 57: 233–38.
    Staniswalis, J. G., K. Messer, and D. R. Finston. 1993. “Kernel Estimators for Multivariate Regression.” Journal of the Nonparametric Statistics, no. 3: 103–21.
    Tapia, R. A. 1971. “The Differentiation and Integration of Nonlinear Operators.” Academic Press, New York.
    Tarter, M. E., and R. A. Kronmal. 1970. “On Multivariate Density Estimates Based on Orthogonal Expansions.” Annals of the Mathematical Statistics, no. 41: 718–22.
    ———. 1976. “An Introduction to the Implementation and Theory of Nonparametric Density Estimation.” American Statistician, no. 30: 105–12.
    Taylor, C. C. 1989. “Bootstrap Choice of the Smoothing Parameter in Kernel Density Estimation.” Biometrika, no. 76: 705–12.
    Terrell, G. R. 1984. “Efficiency of Nonparametric Density Estimators.” Technical Report, Department of Math Sciences, Rice University.
    ———. 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. 1980. “On Improving Convergence Rates for Nonnegative Kernel Density Estimators.” Annals of Statistics, no. 8: 1160–63.
    Terrell, G. R., and D. W. Scott. 1985. “Oversmoothed Nonparametric Density Estimates.” Journal of the American Statistical Association, no. 80: 209–14.
    Terrell, G. R., and D. W. Scott. 1992a. “Variable Kernel Density Estimation.” Annals of Statistics, no. 20: 1236–65.
    ———. 1992b. “Variable Kernel Density Estimation.” Annals of Statistics, no. 20: 1236–65.
    Tierney, L. 1990. LISP-STAT. John Wiley, New York.
    Wahba, G. 1971. “A Polynomial Algorithm for Density Estimation.” Annals of the Mathematical Statistics, no. 42: 1870–86.
    ———. 1977. “Optimal Smoothing of Density Estimates.” Academic Press, New York.
    ———. 1981. “Data-Based Optimal Smoothing of Orthogonal Series Density Estimates.” Annals of Statistics, no. 9: 146–56.
    Walter, G., and J. R. Blum. 1979a. “Probability Density Estimation Using Delta Sequences.” Annals of Statistics, no. 7: 328–40.
    ———. 1979b. “Probability Density Estimation Using Delta Sequences.” Annals of Statistics, no. 7: 328–40.
    Wand, M. P. 1994. “Fast Computation of Multivariate Kernel Estimators.” Journal of Computational and Graphical Statistics, no. 3: 433–45.
    Wand, M. P., and M. C. Jones. 1991a. “Comparison of Smoothing Parameterizations in Bivariate Density Estimation.” Journal of the American Statistical Association, no. 88: 520–28.
    ———. 1991b. “Comparison of Smoothing Parameterizations in Bivariate Density Estimation.” Journal of the American Statistical Association, no. 88: 520–28.
    ———. 1995. Kernel Smoothing. CRC Press, London.
    Watson, G. S. 1969. “Density Estimation by Orthogonal Series.” Annals of the Mathematical Statistics, no. 40: 1496–98.
    Watson, G. S., and M. R. Leadbetter. 1963. “On the Estimation of the Probability Density i.” Annals of the Mathematical Statistics, no. 34: 480–91.
    Woodroofe, M. 1970. “On Choosing a Delta-Sequence.” Annals of the Mathematical Statistics, no. 41: 1665–71.
    Worton, B. J. 1989. “Optimal Smoothing Parameters for Multivariate Fixed and Adaptive Kernel Methods.” Jouranl of the Statistical Computations and Simulations, no. 32: 45–57.