Capítulo 2

Critérios de estimação não paramétricos


O foco da estimação não paramétrica difere do da estimação paramétrica. No caso da estimação paramétrica, dada uma família de densidade paramétrica \(f(\cdot\, | \, \theta )\), como a família normal de dois parâmetros \(N(\mu,\sigma^2)\), onde \(\theta = (\mu,\sigma^2)\), a ênfase recai na obtenção do melhor estimador \(\widehat{\theta}\) de \(\theta\). No caso não paramétrico, a ênfase recai diretamente na obtenção de uma boa estimativa \(\widehat{f}(\cdot)\) da função de densidade completa \(f(\cdot)\). Neste capítulo, apresentamos uma introdução aos critérios de estimação não paramétrica, utilizando apenas ferramentas já conhecidas da análise paramétrica.

R. A. Fisher e Karl Pearson travaram um debate acirrado sobre aspectos da estimação paramétrica, com Pearson defendendo o uso de curvas não paramétricas. As curvas não paramétricas são determinadas pela estrutura dos dados e têm ampla aplicabilidade. As curvas paramétricas dependem da construção de modelos e do conhecimento prévio das equações subjacentes aos dados. Fisher (1922, 1932) denominou as duas fases da estimação paramétrica como os problemas de especificação e estimação. Fisher concentrou-se em alcançar a otimalidade na fase de estimação. Sua única observação sobre o problema da especificação foi que cabia “inteiramente ao estatístico prático” escolher uma forma que “sabemos como lidar” com base na experiência.

Fisher observou que a especificação incorreta poderia ser detectada por um teste a posteriori, mas não ofereceu instruções adicionais. Se um modelo paramétrico for superparametrizado na tentativa de fornecer maior generalidade, então muitas escolhas do vetor \(\theta\) fornecerão estimativas pontuais quase idênticas \(\widehat{f}(x)\). Um modelo paramétrico especificado incorretamente apresenta um viés que não pode ser removido apenas por grandes amostras.

Determinar se o viés é grande demais para justificar a manutenção do modelo paramétrico pode ser complexo, visto que os testes de ajuste quase sempre rejeitam modelos bastante razoáveis com amostras grandes. A “maldição da otimalidade” reside no fato de que a aplicação incorreta de métodos “ótimos” é preferida a métodos mais gerais, porém menos eficientes.

O que Pearson não conseguiu argumentar de forma convincente é que estimadores ótimos podem se tornar ineficientes com pequenas perturbações nas premissas subjacentes ao modelo paramétrico. A ênfase moderna na estimação robusta sacrifica, corretamente, uma pequena porcentagem da otimalidade paramétrica para alcançar maior insensibilidade à especificação incorreta do modelo. Contudo, em muitas situações multivariadas, apenas informações prévias vagas sobre um modelo apropriado estão disponíveis. Os métodos não paramétricos eliminam a necessidade de especificação do modelo. A perda de eficiência não precisa ser muito grande e é compensada pela redução do risco de interpretação errônea dos dados devido à especificação incorreta do modelo.



2.1 Estimação da função de distribuição


A função mais simples de estimar não parametricamente é a função de distribuição acumulada (FDA) de uma variável aleatória \(X\), definida por \[ F(x)=P(X\leq x)\cdot \]

O estimador óbvio da teoria elementar da probabilidade é a função de distribuição acumulada empírica (FDC), definida como \[ \tag{2.1} F_n(x)=\dfrac{\# \{x_i\leq x\}}{n} = \dfrac{\# x_i\in (-\infty,x]}{n}=\dfrac{1}{n}\sum_{i=1}^n \pmb{I}_{(-\infty,x]}(x_i), \] onde \(\{x_1,x_2,\cdots,x_n\}\) é uma amostra aleatória de \(F\) e \[ \pmb{I}_A(x)=\left\{\begin{array}{rcl} 1, & \mbox{se} & x\in A \\[0.em] 0, & \mbox{se} & x\notin A \end{array}\right.\cdot \]

Essa função tem um formato de escada, como mostrado na Figura 2.1 nos dados do gêiser apresentados a continuação.

geyser = c(4.37, 3.87, 4, 4.03, 3.5, 4.08, 2.25, 4.7, 1.73, 4.93, 1.73, 4.62, 3.43,
    NA, 4.25, 1.68, 3.92, 3.68, 3.1, 4.03, 1.77, 4.08, 1.75, 3.2, 1.85, 4.62,
    1.97, NA, 4.5, 3.92, 4.35, 2.33, 3.83, 1.88, 4.6, 1.8, 4.73, 1.77, 4.57,
    1.85, 3.52, NA, 4, 3.7, 3.72, 4.25, 3.58, 3.8, 3.77, 3.75, 2.5, 4.5, 4.1,
    3.7, 3.8, 3.43, NA, 4, 2.27, 4.4, 4.05, 4.25, 3.33, 2, 4.33, 2.93, 4.58,
    1.9, 3.58, 3.73, 3.73, NA, 1.82, 4.63, 3.5, 4, 3.67, 1.67, 4.6, 1.67, 4,
    1.8, 4.42, 1.9, 4.63, 2.93, NA, 3.5, 1.97, 4.28, 1.83, 4.13, 1.83, 4.65,
    4.2, 3.93, 4.33, 1.83, 4.53, 2.03, NA, 4.18, 4.43, 4.07, 4.13, 3.95, 4.1,
    2.72, 4.58, 1.9, 4.5, 1.95, 4.83, 4.12)

Os dados em geyser são a duração em minutos de 107 erupções quase consecutivas do gêiser Old Faithful. Fonte: Weisberg (1985).

library(ggplot2)
sample_geyser <- data.frame(x = geyser)
ggplot(sample_geyser, aes(x)) +
  stat_ecdf(geom = "step") +
  labs(title = "", x = "Duração da erupção (min)", y = expression(paste(F[n],"(x)")))

Figura 2.1: Função de distribuição acumulada empírica do conjunto de dados do gêiser Old Faithful.

Pode-se observar que \(F_n(x)\) possui excelentes propriedades matemáticas para estimar o nível da função \(F(x)\) para cada valor fixo da ordenada \(x\): \[ \mbox{E}\big(F_n(x)\big)=\mbox{E}\big(\pmb{I}_{(-\infty,x]}(X) \big)=1\times P\big(X\in (-\infty,x]\big)=F(x)\cdot \]

De fato, \(nF_n(x)\) é uma variável aleatória binomial \(B(n,p)\), com \(p = F(x)\), de modo que \[ \mbox{Var}\big(F_n(x)\big) = \dfrac{p(1-p)}{n}\cdot. \] Não existem outros estimadores não viesados com variância menor.

Este resultado decorre do fato de que as estatísticas de ordem formam uma estatística suficiente completa, e \(F_n(x)\) é tanto não viesada quanto uma função da estatística suficiente. Observe, porém, que embora a função de distribuição seja frequentemente conhecida por ser contínua, o estimador ótimo \(F_n(x)\) não o é.

Uma pergunta pertinente é se a função de distribuição ou sua função de densidade de probabilidade (fdp) associada, \(f(x) = F'(x)\), deve ser o foco da análise de dados. Certamente, a função de distribuição é bastante útil.

Por um lado, características comuns como assimetria e multimodalidade são muito mais facilmente percebidas em um gráfico da função de densidade do que em um gráfico da função de distribuição, compare as Figuras 1.13 e 2.1. Por outro lado, é difícil ignorar o fato de que muitos cientistas sociais estão mais focados em probabilidades e, portanto, preferem estimar a função de distribuição acumulada (FDA) em vez da função de densidade.

A distinção entre a utilidade da função de distribuição acumulada (CDF) e da função de densidade de probabilidade (fdp) torna-se mais clara em mais de uma dimensão. A definição de uma função de distribuição acumulada empírica multivariada é simplesmente \[ F_n(\pmb{x})=\dfrac{\#\{\pmb{x}_i\leq \pmb{x}\}}{n}, \qquad \pmb{x}\in\mathbb{R}^d, \] onde a desigualdade \(\{\pmb{x}_i\leq \pmb{x}\}\) é aplicada componente a componente para cada ponto de dados. A mesma propriedade de otimalidade se mantém como no caso univariado. No entanto, poucos estatísticos já viram uma função de distribuição cumulativa empírica bivariada.

Considere a função de distribuição acumulada empírica bivariada (ECDF) na Figura 2.2 do conjunto de dados do gêiser, correspondente à função de densidade de probabilidade bivariada (PDF) mostrada na Figura 1.13. Dado o pequeno número de lacunas na série temporal de 107 observações, existem 106 pares \(\{x_t,x_{t+1}\}\) incluídos nesta figura.

library(robusTest)
N = length(geyser)
x = na.omit(geyser[1:(N-1)])
y = na.omit(geyser[2:N])
library(Emcdf)
data = cbind(x,y)
plotcdf(data, type = "wireframe", main = "", angle = 40)

Figura 2.2: Função de distribuição acumulada bivariada empírica dos dados defasados de Old Faithful \(\{x_t,x_{t+1}\}\).

A característica trimodal pode ser reconhecida com um pouco de reflexão, mas não tão facilmente quanto na Figura 1.13. A superfície é talvez inesperadamente complexa, dado o pequeno tamanho da amostra. Em particular, o número de saltos na função é considerado no Exercício 2.1. Assim, a função de distribuição multivariada tem pouco interesse para fins gráficos ou de análise de dados. Além disso, aplicações estatísticas multivariadas ubíquas, como regressão e classificação, dependem da manipulação direta da função de densidade e não da função de distribuição.


2.2 Estimação não paramétrica direta da densidade


Seguindo a relação teórica \(f(x) = F'(x)\), a função de densidade de probabilidade empírica (epdf) é definida como a derivada da função de distribuição acumulada empírica: \[ \tag{2.2} f_n(x)=\dfrac{\mbox{d}}{\mbox{d}x} F_n(x)=\dfrac{1}{n}\sum_{i=1}^n \delta(x-x_i), \] onde \(\delta(t)\) é a função delta de Dirac.

A função de densidade de probabilidade empírica (epdf) é sempre uma densidade uniforme discreta sobre os dados, ou seja, a massa de probabilidade é \(n^{-1}\) em cada ponto de dados. A epdf é como um diagrama de dispersão unidimensional (gráfico de pontos) e é uma estimativa inútil se a densidade for contínua (ver Figura 2.3). A epdf é claramente inferior ao histograma do ponto de vista gráfico. O principal uso moderno da epdf tem sido como densidade de amostragem para o bootstrap (Efron 1982).

Figura 2.3: Um histograma e uma função de densidade de probabilidade empírica (apontando para baixo) do conjunto de dados de espessura da moeda de um centavo americana.

Existe um estimador não viesado de variância mínima uniforme de \(f(x)\)? No primeiro tratamento teórico do assunto, Rosenblatt (1956) provou que tal estimador não existe.

Seja \(X_n\) uma amostra aleatória de tamanho \(n\) de \(f\). Suponha que exista um estimador \(T_n(x; X_n)\) tal que \(\mbox{E}\big(T_n(x; X_n)\big) = f(x)\) para toda \(f\) contínua e para todos os valores de \(x\) e \(n\). Os dados devem aparecer simetricamente em \(T_n(x; X_n)\), pois a estimativa não paramétrica não pode variar com a ordem em que a amostra foi coletada. Agora, para todos os intervalos \((a, b)\), \[ \mbox{E}\left( \int_a^b T_n(x)\mbox{d}x\right) = \int_a^b f(x)\mbox{d}x=F(b)-F(a)=\mbox{E}\big( F_n(b)-F_n(a)\big), \] utilizando o teorema de Fubini para justificar a inversão da ordem dos operadores de esperança e integral.

Tanto \(T_n(x; X_n)\) quanto \(F_n(b)-F_n(a)\) são funções das estatísticas suficientes completas e, como \(F_n(b)-F_n(a)\) é o único estimador simétrico não viesado de \(F(b)-F(a)\), segue-se que \[ F_n(b)-F_n(a)=\int_a^b T_n(x)\mbox{d}x, \] para quase todas as amostras \(\{X_n\}\).

Este resultado é uma contradição, visto que o lado direito da equação é absolutamente contínuo, enquanto o lado esquerdo não o é. Na época, esse resultado foi surpreendente e decepcionante em uma profissão que se acostumara à busca por estimadores ótimos não viesados.

Hoje, estimadores viesados são frequentemente preferidos em situações onde estimadores não viesados estão disponíveis, por exemplo, com estimadores de encolhimento como a regressão ridge (Hoerl and Kennard 1970) e os estimadores de Stein (1956), ou com estimadores de penalização como o LASSO (Tibshirani 1996).


2.3 Critérios de erro para estimadores da densidade


O desejo de comparar diferentes estimadores e tentar identificar o melhor pressupõe a especificação de um critério que possa ser otimizado. A otimalidade não é um conceito absoluto, mas está intimamente ligado à escolha de um critério. A preferência por um critério é em grande parte subjetiva, embora certos argumentos teóricos ou intuitivos possam ser apresentados.

No entanto, a eliminação total do elemento subjetivo da estimação não paramétrica parece indesejável; por exemplo, será demonstrado que a quantidade de ruído em histogramas “ótimos” pode evocar uma resposta negativa com conjuntos de dados muito grandes. No mundo paramétrico, um estimador ótimo provavelmente será ótimo para qualquer propósito relacionado. No mundo não paramétrico, um estimador pode ser ótimo para um propósito e notavelmente subótimo para outro. Esse trabalho extra é um preço a ser pago por trabalhar com uma classe mais geral de estimadores.

Ao aproximar parâmetros com estimadores enviesados, o critério de variância é frequentemente substituído pelo erro quadrático médio (EQM), que é a soma da variância e do quadrado do viés. Para a estimação pontual de uma função de densidade pelo estimador \(\widehat{f}(x)\), \[ \mbox{MSE}\big(\widehat{f}(x)\big)=\mbox{E}\big(\widehat{f}(x)-f(x) \big)^2 = \mbox{Var}\big(\widehat{f}(x)\big)+\mbox{Viés}^2\big(\widehat{f}(x)\big), \] onde \(\mbox{Viés}\big(\widehat{f}(x)\big)=\mbox{E}\big(\widehat{f}(x)\big)-f(x)\).

Esta equação trata o problema de estimação da densidade não paramétrica como um problema padrão de estimação pontual com parâmetro desconhecido \(\theta = f(x)\). Embora tais análises pontuais possam ser interessantes, a ênfase nossa será claramente na estimação e exibição de toda a superfície de densidade.

Para alguns, o critério global mais intuitivamente atraente é a norma \(L_\infty\): \[ \sup_x \left| \widehat{f}(x)-f(x)\right|\cdot \]

No outro extremo do espectro está a norma \(L_1\): \[ \int \left| \widehat{f}(x)-f(x)\right| \mbox{d}x\cdot \] Nenhum desses critérios é tão facilmente manipulável quanto a norma \(L_2\), que neste contexto é denominada erro quadrático integrado (\(\mbox{ISE}\)): \[ \tag{2.3} \mbox{ISE} = \int \left( \widehat{f}(x)-f(x)\right)^2 \mbox{d}x\cdot \] Ainda assim, o erro quadrático integrado é uma variável aleatória complexa que depende da verdadeira função de densidade desconhecida, do estimador específico e do tamanho da amostra; além disso, mesmo com essas três quantidades fixas, o \(\mbox{ISE}\) é uma função da realização específica de \(n\) pontos.

Para a maioria dos propósitos, será suficiente examinar a média do \(\mbox{ISE}\) nessas realizações; ou seja, a média da variável aleatória \(\mbox{ISE}\) ou erro quadrático integrado médio (\(\mbox{MISE}\)): \[ \mbox{MISE} = \mbox{E}\big( \mbox{ISE}\big) = \mbox{E}\left( \int \left( \widehat{f}(x)-f(x)\right)^2 \mbox{d}x\right) = \int \mbox{E}\left( \widehat{f}(x)-f(x)\right)^2 \mbox{d}x = \\[0.8em] =\int \mbox{MISE}\big( \widehat{f}(x)\big)\mbox{d}x = \mbox{IMSE}, \] onde a troca dos operadores integral e esperança é justificada por uma aplicação do teorema de Fubini.

A última grandeza é o \(\mbox{IMSE}\), abreviação de erro quadrático médio integrado. Assim, o critério de erro \(\mbox{MISE}\) possui duas interpretações diferentes, embora equivalentes: ele mede tanto o erro médio global quanto o erro pontual acumulado.

Outros candidatos possíveis para medir o erro incluem números de informação, como o critério de Kullback-Leibler, que é definido como \[ \int f(x)\log\big(f(x)/\widehat{f}(x) \big)\mbox{d}x, \] Distância de Hellinger, que é definida como \[ \left(\int \Big(\widehat{f}^{1/p}(x) - f^{1/p}(x)\Big)^p \mbox{d}x\right)^{1/p}, \] geralmente com \(p = 2\), critério de informação de Akaike, ou mesmo outras distâncias \(L_p\). Essas alternativas não serão exploradas aqui.


2.3.1 Erro quadrático médio integrado (\(\mbox{MISE}\)) para estimadores paramétricos


Poucos estatísticos têm muita intuição sobre o critério \(\mbox{MISE}\), visto que ele quase nunca é discutido no contexto paramétrico. No entanto, é perfeitamente possível avaliar a qualidade dos estimadores de densidade paramétricos pelo \(\mbox{MISE}\) da estimativa de densidade paramétrica completa, em vez do \(\mbox{MSE}\) do parâmetro isoladamente. Na maioria dos casos, tanto o \(\mbox{MISE}\) do parâmetro quanto o \(\mbox{MSE}\) da estimativa de densidade diminuem a uma taxa de \(O(n^{-1})\) à medida que o tamanho da amostra aumenta.


Exemplo 2.1: Densidade uniforme. considere a função de densidade uniforme \(f=U(0,\theta)\), onde \(\theta\) vai ser estimado pelo estimador de máxima verossimilhança \(\widehat{\theta}=x_{(n)}\), a \(n\)-ésima estatística de ordem. Assim \[ \widehat{f}=U(0,x_{(n)})\cdot \] Sem perda de generalidade, escolhemos \(\theta=1\). Seguindo a equação (2.3), \[ \mbox{ISE}=\left(\dfrac{1}{x_{(n)}}-1 \right)^2 x_{(n)}+(0-1)^2 (1-x_{(n)})=\dfrac{1}{x_{(n)}}-1\cdot \] Dado que, \[ f(x_{(n)})=nx_{(n)}^{n-1}, \] para \(0<x_{(n)}<1\), segue que \[ \tag{2.4} \mbox{MISE}=\int_0^1 \left(\dfrac{1}{x_{(n)}}-1 \right)^2 nx_{(n)}^{n-1}\mbox{d}x_{(n)}=\dfrac{1}{n-1}=O(n^{-1})\cdot \]

Considere a família unidimensional de estimadores \(\widehat{f}=U(0,cx_{(n)})\) , onde \(c\) é uma constante positiva. Existem dois casos, dependendo se \(cx_{(n)} < 1\) ou não, mas o \(\mbox{ISE}\) pode ser calculado em cada instância (veja o Exercício 2.2).

Tomando as esperanças, obtemos: \[ \tag{2.5} \mbox{MISE}(c)=\left\{ \begin{array}{ccl} \dfrac{n}{(n-1)c}-1, & \mbox{caso} & c<1 \\[0.8em] \dfrac{2-nc^{n-1}+(n-1)c^n}{(n-1)c^n}, & \mbox{caso} & c>1 \end{array}\right.\cdot \]

Como ambos os \(\mbox{MISE}(c)\) e suas derivadas são contínuas em \(c = 1\), com \[ \mbox{MISE}'(c = 1) = -\dfrac{n}{(n-1)}, \] o \(\mbox{MISE}\) mínimo é obtido para algum \(c > 1\), ou seja, \[ c^*=2^{1/(n-1)}\approx 1+n^{-1}\log(2)\cdot \] Este resultado é semelhante ao fato bem conhecido de que o \(\mbox{MSE}\) do estimador paramétrico é minimizado quando \(c = 1 + (n+1)^{-1}\). No entanto, neste caso, o \(\mbox{MSE}\) paramétrico converge a uma taxa excepcionalmente rápida de \(O(n^{-2})\) (ver Romano and Siegel 1986).

A taxa do \(\mbox{MISE}\) de \(O(n^{-1})\) reflete com mais precisão o erro médio para toda a curva de densidade e não apenas para seu ponto final.


Exemplo 2.2: Método \(\mbox{MISE}\) paramétrico geral com aplicação gaussiana. Considere um estimador não viesado \(\widehat{\theta}\) de um parâmetro \(\theta\). Escrevendo o \(\mbox{ISE}\) paramétrico como \(\pmb{I}(\widehat{\theta})\)e tomando uma série de Taylor, obtemos \[ \pmb{I}(\widehat{\theta})=\int \Big(f(t\, | \, \widehat{\theta})-f(t\, | \, \theta) \Big)^2 \mbox{d}t =\sum_k \dfrac{1}{k!} \big(\widehat{\theta}-\theta \big)^k \pmb{I}^{(k)} (\theta)\cdot \] Agora \(\pmb{I}(\theta)=0\) e \[ \pmb{I}'(\theta)=2 \int \Big(f(t\, | \, \widehat{\theta})-f(t\, | \, \theta) \Big) f'(t\ | \, \theta) \mbox{d}t = 0 \] também; portanto, \[ \mbox{MISE}(\theta) = \mbox{E}\big( \pmb{I}(\widehat{\theta})\big)=\dfrac{1}{2}\mbox{Var}\big(\widehat{\theta}\big)\pmb{I}''(\widehat{\theta})+\cdots \; \cdot \] A omissão de termos de ordem superior na expansão do \(\mbox{MISE}\) resulta no erro quadrático médio integrado assintótico (\(\mbox{AMISE}\)), que é dado por \[ \mbox{AMISE}(\theta)=\dfrac{1}{2}\pmb{I}''(\widehat{\theta})\mbox{Var}\big(\widehat{\theta}\big) \cdot \] Este resultado pode ser estendido a um vetor de parâmetros por meio de uma expansão de Taylor multivariada.

Considere a estimação da densidade normal de dois parâmetros. \[ \phi(x\, | \, \mu,\sigma^2)=\dfrac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\dfrac{(x-\mu)^2}{2\sigma^2} \right), \] e um estimador da forma \[ \widehat{\phi}(x)=\phi(x\, |\, u,\nu); \] Por exemplo, \(\widehat{\phi}(x) = \phi(x\, | \, \overline{x},s^2)\), usando estimativas não viesadas de \(\mu\) e \(\sigma^2\). Se a densidade verdadeira for uma normal padrão, então o \(\mbox{ISE}\) bivariado pode ser expandido como \[ \tag{2.6} \pmb{I}(u,\nu)\approx \pmb{I}(0,1)+u\pmb{I}_u+(\nu-1)\pmb{I}_\nu +\dfrac{u^2}{2}\pmb{I}_{uu}+u(\nu-1)\pmb{I}_{u\nu}+\dfrac{(\nu-1)^2}{2}\pmb{I}_{\nu\nu}, \] onde as derivadas parciais \((\pmb{I}_u,\pmb{I}_\nu,\pmb{I}_{uu},\pmb{I}_{u\nu},\pmb{I}_{\nu\nu})\) são avaliadas em \((u,\nu)=(0,1)\).

Então (veja o Exercício 2.4) \[ \tag{2.7} \begin{array}{rcl} \pmb{I}_u(u,\nu) & = & \displaystyle 2\int (\widehat{\phi}-\phi)\left( \dfrac{x-u}{\nu}\widehat{\phi}\right)\mbox{d}x, \\[0.8em] \pmb{I}_{uu}(u,\nu) & = & \displaystyle 2\int (\widehat{\phi}-\phi)\left( \dfrac{x-u}{\nu}\widehat{\phi}\right)'\mbox{d}x + 2\int \left( \dfrac{x-u}{\nu}\widehat{\phi}\right)^2\mbox{d}x \cdot \end{array} \]

Agora, o termo \((\widehat{\phi}-\phi)\) quando avaliado em \((u,\nu)=(0,1)\). A integral final torna-se \[ 2\int \big( x\phi(x\, | \, 0,1)\big)^2\mbox{d}x =\dfrac{1}{2\sqrt{\pi}}\cdot \] De maneira similar \(\pmb{I}_{\nu\nu}(0,1)=3/(16\sqrt{\pi})\) e \(\pmb{I}_{u\nu}(0,1)=0\). Portanto, tomando a esperança da equação (2.6), obtemos \[ \tag{2.8} \mbox{AMISE}=\dfrac{1}{2}\pmb{I}_{uu}\mbox{Var}(u)+\dfrac{1}{2}\pmb{I}_{\nu\nu}\mbox{Var}(\nu) \approx \dfrac{1}{4n\sqrt{\pi}}+\dfrac{3}{16n\sqrt{\pi}}=\dfrac{7}{16n\sqrt{\pi}}, \] visto que as variações de \(\overline{x}\) e \(s^2\) são \(1/n\) e \(2/(n-1)\approx 2/n\), respectivamente. É interessante notar que o \(\mbox{AMISE}\) na equação (2.8) é a soma de \[ \tag{2.9} \mbox{AMISE}\big(\phi(x\, | \, \overline{x},1)\big)+\mbox{AMISE}\big(\phi(x\, | \, 0,s^2)\big)=\dfrac{1}{4n\sqrt{\pi}}+\dfrac{3}{16n\sqrt{\pi}} \approx \dfrac{0.2468}{n}\cdot \]

Portanto, o \(\mbox{MISE}\) paramétrico é maior se a média for desconhecida do que se a variância for desconhecida, porque as curvas normais deslocadas estão mais distantes umas das outras do que quando as curvas normais se cruzam, como quando apenas suas variâncias diferem. Este cálculo simples foi verificado por uma estimativa de simulação de Monte Carlo de \(1/(4n)\), que é apenas 1.2% maior.


2.3.2 O critério \(L_1\)


2.3.2.1 \(L_1\) vesus \(L_2\)


Intuitivamente, um dos atrativos do critério \(L_1\) \(\int |\widehat{f}-f|\mbox{d}x\) é que ele dá mais atenção às caudas de uma função de densidade do que o critério \(L_2\), que minimiza os valores de densidade relativamente pequenos ao elevá-los ao quadrado. Essa minimização pode ser verificada em casos específicos e por simulação. Teoricamente, o critério \(L_1\) apresenta diversas outras vantagens importantes. Primeiro, considere uma análise dimensional.

Como uma função de densidade tem o inverso do comprimento como unidade, \(L_1\) é uma grandeza adimensional após a integração. O critério \(L_2\), por outro lado, mantém as unidades do inverso do comprimento após a integração do erro quadrático da densidade. É possível tentar tornar \(L_2\) adimensional de maneira semelhante à construção de tais grandezas para a assimetria e curtose, mas nenhum método totalmente satisfatório está disponível (ver Seção 7.2.1).

Além disso, \(L_1\) é invariante a mudanças contínuas e monotônicas de escala, como mostrado abaixo. Também pode-se ver que \(0\leq L_1\leq 2\), enquanto \(0\leq L_2\leq \infty\) (ver Exercício 2.7). Por outro lado, em situações práticas, os estimadores que otimizam esses critérios são semelhantes. O ponto a ser lembrado é que não existe um estimador canônico para o método \(L_1\) ou para o método \(L_2\).

Devroye and Györfi (1985) compilaram um tratamento teórico impressionante baseado no erro \(L_1\), mas a simplicidade analítica do erro quadrático e sua adequação em aplicações práticas o tornam o critério de escolha aqui. Alguns resultados assintóticos para estimativas \(L_1\), obtidos por Hall and Wand (1988a) e Scott and Wand (1991), corroboram a noção de que as diferenças práticas entre os critérios \(L_1\) e \(L_2\) são razoavelmente pequenas, exceto em situações extremas.


2.3.2.2 Três propriedades úteis do critério \(L_1\)


Vários aspectos teoricamente interessantes do erro absoluto são discutidos em detalhes mais adiante. Cada um deles também se aplica ao contexto multivariado (Devroye 1987).

A primeira e mais atraente propriedade é a sua interpretabilidade, visto que \(L_1\) é adimensional e invariante sob qualquer transformação monótona suave. Portanto, é possível comparar a dificuldade relativa de estimar diferentes densidades.

Para ver isso, suponha \(X\sim f\) e \(Y\sim g\), definido \(X^*=h(X)\) e \(Y^*=h(Y)\) temos \[ f^* = f\big(h^{-1}(x^*) \big)|J| \] com uma expressão similar para \(g^*\), onde \(J\) é o Jacobiano.

Uma simples mudança de variáveis resulta em \[ \begin{array}{rcl} \displaystyle \int_u |f^*(u)-g^*(u) |\mbox{d}u & = & \displaystyle \int_u \left| f\big(h^{-1}(u) \big)-g\big(h^{-1}(u) \big)\right| |J|\mbox{d}u\\[0.8em] & = & \displaystyle \int_\nu \left| f(\nu)-g(\nu)\right| \mbox{d}\nu\cdot \end{array} \] Uma interpretação incorreta desse resultado é concluir que todas as escalas de medição são igualmente difíceis de estimar ou que a transformação de variáveis não afeta a qualidade de uma estimativa de densidade.

O segundo resultado é conhecido como Lema de Scheffé: \[ \tag{2.10} \begin{array}{rcl} 2\sup_A \left|\int_A f(x)\mbox{d}x -\int_A g(x)\mbox{d}x\right| & = & \displaystyle \int_A \left| f(x)\mbox{d}x - g(x)\right|\mbox{d}x \\[0.8em] & = & \displaystyle 2\int _{\{x\in A\, : \, f(x)>g(x)\}} \big(f(x)-g(x)\big)\mbox{d}x\cdot \end{array} \]

Para provar isso, considere o conjunto \(B=\{x\,: \, f(x)>g(x)\}\); então para qualquer conjunto \(A\), \[ \begin{array}{rcl} \displaystyle 2 \int_A \big( f(x) - g(x)\big) \mbox{d}x & = & \displaystyle 2\int_{A\cap B} \big( f(x) - g(x)\big) \mbox{d}x -2 \int_{A\cap B^c} \big( g(x) - f(x)\big) \mbox{d}x \leq 2\int_B \big( f(x) - g(x)\big) \mbox{d}x - 0 \\[0.8em] & = & \displaystyle \int_B \big( f(x) - g(x)\big) \mbox{d}x + \left(1-\int_{B^c} f(x)\mbox{d}x \right) - \left( 1-\int_{B^c} g(x)\mbox{d}x \right) \\[0.8em] & = & \displaystyle \int_B \big( f(x) - g(x)\big) \mbox{d}x + \int_{B^c} \big( g(x) - f(x)\big) \mbox{d}x = \int \left| f(x) - g(x) \right| \mbox{d}x, \end{array} \] o que estabelece a segunda igualdade no lema; tomando o supremo sobre \(A\), estabelece-se o “\(\leq\)” para a primeira igualdade. O “\(\geq\)” segue diretamente de \[ 2\sup_A \left| \int_A \big( f(x) - g(x)\big) \mbox{d}x \right|\geq 2\int_{A=B} \big( f(x) - g(x)\big) \mbox{d}x\cdot \]

A igualdade (2.10) estabelece uma conexão com a classificação estatística. Suponha que os dados provenham aleatoriamente de duas densidades, \(f\) e \(g\), e que um novo ponto \(x\) deva ser classificado. Usando uma regra Bayesiana da forma: atribua \(x\) a \(f\) se \(x\in A\), para algum conjunto \(A\), e a \(g\) caso contrário, a probabilidade de classificação incorreta é \[ \begin{array}{rcl} P(\mbox{erro}) & = & \displaystyle \dfrac{1}{2}\int_{A^c} f(x)\mbox{d}x+\dfrac{1}{2}\int_A g(x)\mbox{d}x = \dfrac{1}{2}\left(1-\int_A f(x)\mbox{d}x \right)+\dfrac{1}{2}\int_A g(x)\mbox{d}x \\[0.8em] & = & \displaystyle \dfrac{1}{2}-\dfrac{1}{2}\int_A \big( f(x) - g(x)\big) \mbox{d}x\cdot \end{array} \]

Escolher \(A\) para minimizar esse erro leva a \(A = B\), usando o lema, e nos dá o terceiro resultado \[ \tag{2.11} P(\mbox{erro}) = \dfrac{1}{2}-\dfrac{1}{4}\int \big( f(x) - g(x)\big) \mbox{d}x\cdot \]

Assim, minimizar a distância \(L_1\) entre \(g=\widehat{f}\) e \(f\) é equivalente a maximizar a probabilidade em (2.11); isto é, otimizar \(L_1\) equivale a maximizar a confusão entre \(\widehat{f}\) e \(f\). Essa otimização é precisamente o que se deseja. A interpretação probabilística é válida em qualquer dimensão. No Exercício 2.8, essa expressão é verificada nos dois casos extremos.


2.3.3 Critérios de estimação paramétrica amostrais


Nesta seção, o estimador paramétrico será denotado por \(f_\theta(x)\), mas a densidade verdadeira será denotada por \(g(x)\). Isso enfatizará que o modelo é apenas uma aproximação da densidade verdadeira. Qual é a base teórica da máxima verossimilhança?

A resposta é a divergência KL (Kullback and Leibler 1951) entre \(f_\theta\) e \(g\): \[ \tag{2.12} \begin{array}{rcl} d_{KL}(f_\theta,g) & = & \displaystyle \int g(x)\log\left( \dfrac{g(x)}{f_\theta(x)}\right)\mbox{d}x\\[0.8em] & = & \displaystyle \int g(x)\log\big(g(x)\big)\mbox{d}x -\int g(x)\log\big(f_\theta(x) \big)\mbox{d}x, \end{array} \] que é não negativa por aplicação da Desigualdade de Jensen (Lehmann and Casella 1998) e é zero se, e somente se, \(f_\theta = g\).

A primeira integral em (2.12) é constante, e tentamos substituir a densidade verdadeira desconhecida \(g(x)\) pela função de densidade de probabilidade empírica \(f_n(x)=\frac{1}{n}\sum_{i=1}^n \delta(x-x_i)\), ver equação 2.2.

Agora \(\int \delta(x-x_i)\log\big(f_\theta(x)\big)\mbox{d}x=\log\big(f_\theta(x_i)\big)\); então \[ \widehat{d}_{KL}(f_\theta,g)=\mbox{constante}-\dfrac{1}{n}\sum_{i=1}^n \log\big(f_\theta(x_i)\big)\cdot \] Assim, minimizar a divergência de Kullback-Leibler (KL) é equivalente a maximizar o logaritmo da verossimilhança.

Em seguida, considere um dos critérios adimensionais, a saber, a Distância de Hellinger: \[ d_H(f_\theta,g)^2 = \int \left(\sqrt{f_\theta(x)}-\sqrt{g(x)} \right)^2 \mbox{d}x\geq 0\cdot \] Não existe uma grandeza óbvia para substituir a incógnita \(g(x)\), exceto uma estimativa não paramétrica, como um histograma.

Além da complexidade computacional da minimização numérica necessária, diferentes escolhas do histograma resultariam em diferentes estimativas de \(\theta\). Em contraste, o uso da função de densidade de probabilidade amostral (epdf) em vez de um histograma na divergência de Kullback-Leibler (KL) fornece um critério bem definido e totalmente baseado nos dados. Essas mesmas observações se aplicam ao critério \(L_1\) adimensional.

Finalmente, consideramos um critério que não é adimensional, a saber, o erro quadrático integrado: \[ \tag{2.13} \begin{array}{rcl} \mbox{ISE}(\theta) & = & \displaystyle \int \left( f_\theta(x)-g(x)\right)^2 \mbox{d}x \\[0.8em] & = & \displaystyle \int f_\theta(x)^2\mbox{d}x - 2\int f_\theta(x)g(x)\mbox{d}x + \int g(x)^2\mbox{d}x\\[0.8em] & = & \displaystyle \int f_\theta(x)^2\mbox{d}x - 2\mbox{E}\big(f_\theta(X) \big)+\mbox{constante}\cdot \end{array} \] Assumindo que o modelo seja de quadrado integrável em uma forma conveniente e escolhendo o estimador não viesado óbvio para o termo de esperança, chegamos ao critério totalmente baseado em dados \[ \tag{2.14} \widehat{\theta}=\arg \min_\theta \widehat{\mbox{ISE}}(\theta) = \arg \min_\theta \left(\int f_\theta(x)^2 \mbox{d}x -\dfrac{2}{n}\sum_{i=1}^n f_\theta(x_i) \right)\cdot \] Este critério de distância mínma \(L_2\) foi denominado \(L_2 E(\theta)\) por Scott (2001) e é um caso especial de uma família de divergência investigada por Basu et al. (1998).

Como exemplo, considere o ajuste de uma densidade normal, \(\phi(x \, | \, \mu,\sigma^2)\), aos dados Rayleigh (\(n = 15\)) que mediram o peso de um volume padrão de nitrogênio (ver Tukey 1977). O critério (2.14) a ser minimizado em \(\theta = (\mu,\sigma)\) é \[ \tag{2.15} \widehat{\theta}=\arg \min_\theta L_2 E(\theta) = \arg \min_\theta \left(\dfrac{1}{2\sqrt{\pi\sigma}}-\dfrac{2}{n}\sum_{i=1}^n \phi(x_i \, | \, \mu,\sigma^2) \right)\cdot \] O quadro à esquerda da Figura 2.4 exibe um histograma desse pequeno conjunto de dados e os dois ajustes normais. Uma análise cuidadosa dos dados ao longo do eixo mostra dois agrupamentos (clusters).

O modelo de mistura de quatro parâmetros, \[ \tag{2.16} f_\theta(x)= \omega \phi(x\, | \, \mu_1,\sigma^2)+ (1-\omega)\phi(x \, | \, \mu_2,\sigma^2), \] foi ajustado pelo método \(L_2 E\) e é mostrado no quadro à direita da Figura 4.2, juntamente com um histograma bimodal diferente.

Figura 2.4: (Quadro da esquerda) Histograma com ajustes de MLE e \(L_2 E\) normal aos dados Rayleigh. (Quadro da direita) Ajuste de mistura \(L_2 E\) normal aos dados Rayleigh borrados com variância comum.

Os dados foram suavizados, pois o método \(L_2 E\) é sensível a empates nos dados. Observe que os dois valores suavizados em torno de 2.3013 não estão incluídos na mistura à esquerda (tratados como outliers). O reconhecimento do segundo cluster por Lord Rayleigh resultou na concessão do Prêmio Nobel de Física de 1904 pela descoberta do gás nobre argônio.

Note que o critério \(\mbox{ISE}/L_2 E\) é intuitivamente satisfatório, pois o objetivo é aproximar ao máximo as duas curvas. Consequentemente, focaremos principalmente nesse critério para nossos estimadores de curvas não paramétricos. Donoho and Liu (1988) mostraram que o critério \(\mbox{ISE}\) é robusto por natureza, uma propriedade compartilhada por todos os critérios de distância mínima. Os dois pontos ignorados pelo método \(L_2 E\) no ajuste da mistura são uma evidência dessa propriedade.


2.4 Famílias de distribuições não paramétricas


2.4.1 Família de distribuições de Pearson


Karl Pearson impulsionou o desenvolvimento de estimadores não paramétricos de duas maneiras. Ele cunhou o termo histograma e estudou as funções de densidade que são soluções da equação diferencial \[ \tag{2.17} \dfrac{\mbox{d}\log\big(f(x) \big)}{\mbox{d}x}=\dfrac{x-a}{b+cx+dx^2}\cdot \] Pearson (1902a), Pearson (1902b) identificou sete tipos de soluções para esta equação, dependendo das raízes do denominador e de quais parâmetros eram 0. Curiosamente, esta classe contém a maioria das distribuições clássicas importantes: Normal, \(t\)-Student, Beta e \(F\) de Snedecor (ver Exercício 2.9). Pearson propôs usar os quatro primeiros momentos amostrais para estimar as constantes desconhecidas \((a, b, c, d)\) na equação (2.17). Hoje, a máxima verossimilhança poderia ser usada.

Pearson motivou a equação diferencial (2.17) recorrendo à distribuição hipergeométrica discreta, que descreve uma urna grande com \(N\) bolas, das quais \(pN\) são pretas. Uma amostra de \(n\) bolas é retirada e o número de bolas pretas é registrado como \(X\); então \[ f(x)=P(X=x)=\binom{Np}{x}\binom{N(1-p)}{n-x}\Big/ \binom{N}{n}, \qquad x=0,\cdots,N\cdot \]

A equação diferencial de Pearson (2.17) surge após alguns cálculos (ver Exercício 2.11), mantendo \(p\) fixo à medida que a urna cresce: \[ \tag{2.18} \dfrac{\mbox{d}\log\big(f(x) \big)}{\mbox{d}x}\approx \dfrac{\Delta f(x)}{f(x)}=\dfrac{f(x)-f(x-1)}{f(x)}=\dfrac{x-a}{b+cx+dx^2}, \] onde \[ \tag{2.19} a=b=\dfrac{-(n+1)(Np+1)}{N+2}; \quad c=\dfrac{Np+n+2}{N+2}; \quad d=\dfrac{-1}{N+2}\cdot \]

A contribuição e a perspicácia de Pearson consistiram em usar os dados não para ajustar uma forma paramétrica específica de densidade, mas sim para calcular os coeficientes da equação diferencial a partir dos dados por meio dos momentos amostrais.

Matematicamente, a família de Pearson poderia ser considerada paramétrica; porém, filosoficamente, é não paramétrica. Pearson dedicou recursos consideráveis ao cálculo de percentis para suas distribuições. Dada a falta de poder computacional da época, é notável que a família de Pearson tenha obtido ampla aceitação. Ela ainda aparece regularmente em aplicações práticas atualmente.

Muitas outras famílias de distribuições, algumas multivariadas, foram propostas. A família de Johnson (1949) é notável, juntamente com as generalizações multivariadas de Marshall and Olkin (1985).


2.4.2 Quando um estimador é não paramétrico?


Um estimador da densidade paramétrico é definido pelo modelo \(\widehat{f}(x\, | \, \theta)\), onde \(\theta\in\Theta\). Tem-se mostrado surpreendentemente difícil formular uma definição funcional para o que constitui um estimador de densidade não paramétrico.

Uma definição heurística pode ser proposta com base na condição necessária de que o estimador “funcione” para uma “grande” classe de densidades verdadeiras. Uma noção útil é que um estimador não paramétrico deve ter muitos parâmetros, na verdade, talvez um número infinito, ou um número que diverge em função do tamanho da amostra.

A família de Pearson não se qualifica completamente como não paramétrica sob nenhuma dessas definições, embora a dimensão da família de Pearson seja maior do que a da maioria das famílias paramétricas. Tapia and Thompson (1978) tendem a favorecer a noção de que o estimador não paramétrico deve ser de dimensão infinita. Silverman (1986) simplesmente indica que uma abordagem não paramétrica faz “suposições menos rígidas… sobre a distribuição dos dados observados”. Mas quantos parâmetros um histograma tem?

E quanto ao estimador ingênuo de séries ortogonais com um número infinito de termos que é igual (em distribuição) à função de densidade empírica para qualquer tamanho de amostra? (Veja a Seção 6.1.3.)

Uma definição surpreendentemente elegante está implícita no trabalho de Terrell (Terrell and Scott 1992), que demonstra que todos os estimadores, paramétricos ou não paramétricos, são estimadores de kernel generalizados, pelo menos assintoticamente (ver Seção 6.4).

Terrell introduz a ideia da influência de um ponto de dados \(x_i\) na estimativa de densidade pontual em \(x\). Se \(\widehat{f}(x)\) é um estimador não paramétrico, a influência de um ponto deve se anular assintoticamente se \(|x-x_i| > \epsilon\) para qualquer \(\epsilon > 0\), enquanto a influência de pontos distantes não se anula para um estimador paramétrico. Grosso modo, os estimadores não paramétricos são assintoticamente locais, enquanto os estimadores paramétricos não o são. Contudo, os estimadores não paramétricos não devem ser muito locais para serem consistentes.


2.5 Exercícios


  • 1- Elabore um algoritmo para mostrar graficamente a função de distribuição acumulada bivariada empírica. Indique um intervalo para o número possível de saltos que podem ocorrer nessa função. Forneça exemplos simples para os dois casos extremos. Qual é a ordem do seu algoritmo, ou seja, o número de operações aritméticas em função do tamanho da amostra?


  • 2- Verifique a equação (2.5). Mostre-a graficamente para vários tamanhos de amostra e compare o minimizador real com a fórmula assintótica.


  • 3- Mostre que a distância de Kullback-Leibler esperada para o estimador paramétrico \(\widehat{f} = U(0,x_{(n)})\) de \(f = U(0, 1)\) é \(1/(n-1)\).


  • 4- Verifique as equações em (2.7) e que \(\pmb{I}_{v\nu} (0,1) = 3/(16\sqrt{\pi})\).


  • 5- Complete os cálculos para o \(\mbox{AMISE}\) paramétrico normal. Encontre os estimadores \(\mbox{AMISE}\) ótimos da forma \(c\overline{x}\) e \(cs^2\) para dados normais padrão.


  • 6- Assumindo dados normais padrão, calcule o \(\mbox{MISE}\) exato do estimador \(N(\overline{x},1)\) e verifique se as aproximações em série coincidem. Problema de pesquisa: Você consegue encontrar expressões \(\mbox{MISE}\) de forma fechada semelhantes para os estimadores \(N(0,s^2)\) e \(N(\overline{x},s^2)\)?


  • 7- Verifique se os intervalos dos erros \(L_1\) e \(L_2\) são \((0,2)\) e \((0,\infty)\), respectivamente. Dê exemplos onde os valores extremos são observados.
  • 8- Verifique se a probabilidade de erro dada em (2.11) está correta nos dois casos extremos.


  • 9- Verifique diretamente se as densidades paramétricas mencionadas no texto realmente satisfazem a equação diferencial para o sistema de Pearson.


  • 10- Use um programa simbólico para verificar se \[ f(x) = (1+x^2)^{-k} \exp\big(-\alpha \arctan(x)\big) \] é uma possível solução para a equação diferencial de Pearson. Trace o gráfico de \(f(x)\).


  • 11- Verifique os cálculos nas equações (2.18) e (2.19).


  • 12- Verifique o critério \(L_2 E\) na equação (2.15). Derive o critério \(L_2 E\) de mistura de quatro parâmetros usado na Figura 2.4. Use a função nlminb do R para encontrar \(\widehat{\theta}\). Dica: Use a identidade \[ \int \phi(x\, | \, \mu_1,\sigma_1^2)\phi(x\, | \, \nu_2,\sigma_2^2 )\mbox{d}x = \phi(0 \, | \, \mu_1-\mu_2,\sigma_1^2+\sigma_2^2)\cdot \]


  • 13- Considere uma mistura normal de cinco parâmetros como a equação (2.16) com \(\sigma^2\) substituído por \(\sigma_1^2\) e \(\sigma_2^2\). Suponha que \(\mu_2 = x_1\). Mostre que a verossimilhança diverge para \(+\infty\) quando \(\sigma_2\to 0\). Derive o critério \(L_2 E\) na equação (2.14) para o modelo de mistura de cinco parâmetros. Qual é o limite do critério \(L_2 E\) para este cenário?

2.6 Bibliografia


Basu, A., I. R. Harris, N. L. Hjort, and M. C. Jones. 1998. “Robust and Efficient Estimation by Minimising a Density Power Divergence.” Biometrika, no. 85: 549–59.
Devroye, L. 1987. A Course in Density Estimation. Birkhäuser, Boston.
Devroye, L., and L. Györfi. 1985. Nonparametric Density Estimation: The L1 View. John Wiley, New York.
Donoho, D. L., and R. C. Liu. 1988. “The ‘Automatic’ Robustness of Minimum Distance Functional.” Annals of Statistics, no. 16: 552–86.
Efron, B. 1982. The Jackknife, Bootstrap, and Other Resampling Plans. SIAM, Philadelphia.
Hall, P., and M. P. Wand. 1988a. “Minimizing L1 Distance in Nonparametric Density Estimation.” Journal of the Multivariate Analysis, no. 26 (1988a): 59–88.
Hoerl, A. E., and R. W. Kennard. 1970. “Ridge Regression: Biased Estimation for Non-Orthogonal Problems.” Technometrics, no. 12: 55–67.
Johnson, N. L. 1949. “Systems of Frequency Curves Generated by Methods of Translation.” Biometrika, no. 36: 149–76.
Kullback, S., and R. A. Leibler. 1951. “On Information and Sufficiency.” Annals of Mathematical Statistics, no. 22: 79–86.
Lehmann, E. L., and G. Casella. 1998. Theory of Point Estimation. Second Edition. Springer, New York.
Marshall, A. W., and I. Olkin. 1985. “A Family of Bivariate Distributions Generated by the Bivariate Bernoulli Distribution.” Journal of the American Statistical Association, no. 80: 332–38.
Pearson, K. 1902a. “On the Systematic Fitting of Curves to Observations and Measurements, i.” Biometrika, no. 1: 265–303.
———. 1902b. “On the Systematic Fitting of Curves to Observations and Measurements, II.” Biometrika, no. 2: 1–23.
Romano, J. P., and A. F. Siegel. 1986. Counterexamples in Probability and Statistics. Wadsworth & Brooks/Cole, Pacific Grove, CA.
Rosenblatt, M. 1956. “Remarks on Some Nonparametric Estimates of a Density Function.” Annals of Mathematical Statistics, no. 27: 832–37.
Scott, D. W. 2001. “Parametric Statistical Modeling by Minimum Integrated Square Error.” Technometrics, no. 43: 274–85.
Scott, D. W., and M. P. Wand. 1991. “Feasibility of Multivariate Density Estimates.” Biometrika, no. 78: 197–206.
Silverman, B. W. 1986. Density Estimation for Statistics and Data Analysis. Chapman; Hall, London.
Stein, C. M. 1956. “Inadmissibility of the Usual Estimator for the Mean of a Multivariate Normal Distribution.” Edited by Vol. 1. University of California Press, Berkeley, CA.
Tapia, R. A., and J. R. Thompson. 1978. Nonparametric Probability Density Estimation. John Hopkins University Press, Baltimore.
Terrell, G. R., and D. W. Scott. 1992. “Variable Kernel Density Estimation.” Annals of Statistics, no. 20: 1236–65.
Tibshirani, R. J. 1996. “Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society, B, no. 58: 267–88.
Tukey, J. W. 1977. Exploratory Data Analysis. Addison-Wesley, Reading, MA.
Weisberg, S. 1985. Applied Linear Regression. John Wiley, New York.