De certa forma, problemas de estimação não-paramétrica são extensões de problemas de estimação paramétrica, mas a natureza do primeiro é bem diferente do último. Considere, por exemplo, a situação de \(n\) observações independentes identicamente distribuídas, digamos \(X_1,X_2,\cdots,X_n\).

Em um problema paramétrico, assumimos que a distribuição de \(X_i\) é \(F(\cdot;\theta)\), a qual é totalmente especificada até o vetor de parâmetros \(\theta\); então o problema é essencialmente a estimação de \(\theta\). Em um problema não-paramétrico, a distribuição é totalmente desconhecida com, talvez, algumas restrições em propriedades gerais e, portanto, é denotada por \(F\).

Aqui consideramos estimadores de \(F\) em termos de função de densidade \(f\). A função de densidade tem a vantagem de fornecer uma representação visualmente mais informativa do que a função de distribuição. Por exemplo, o histograma geralmente dá uma ideia aproximada da forma da distribuição. Este último ficou como o único estimador de densidade não paramétrico até 1950. Por essa razão, nossa discussão começará com os histogramas.



1. Histograma


Embora o histograma seja usado extensivamente, não é tão frequente que seja necessária uma definição matemática. Uma maneira de defini-lo é através da função de densidade empírica.


Definição 1: Função de densidade empírica.


Seja \(f\) a derivada de \(F\), por isso, pode-se escrever como \[ f(x) = \lim_{h\to 0} \dfrac{F(x+h)-F(x-h)}{2h}\cdot \]

Dizemos então que \(\widehat{f}\), defnido por \[ \widehat{f}_n(x)=\dfrac{\widehat{F}_n(x+h)-\widehat{F}_n(x-h)}{2h}, \] é o histograma, sendo que \(\widehat{F}\)_n é a função de distribuição empírica.


Um histograma é uma representação gráfica da função de probabilidades ou da função de densidade de um conjunto de dados independentes e foi introduzido pela primeira vez por Karl Pearson no artigo Pearson, K. (1895). Contributions to the Mathematical Theory of Evolution. II. Skew Variation in Homogeneous Material. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 186: 343-414.

A representação mais comum do histograma é um gráfico de barras verticais. A palavra histograma é de origem grega, derivada de duas: histos que pode significar testemunha no sentido de aquilo que se vê, como as barras verticais do histograma e da também palavra grega gramma que significa desenhar, registrar ou escrever.

Para construir um exemplo controlado do gráfico de histograma, simulamos uma amostra de tamanho 150 da distribuição normal padrão, com o comando rnorm e depois construímosum gráfico colorido, abaixo a esquerda. Posteriormente, acrescentamos a este gráfico uma linha com a densidade normal, abaixo a direita.

set.seed(8476)
exemplo = rnorm(150)
par(mar=c(5,4,2,1), mfrow=c(1,2))
hist(exemplo, breaks=12, col="red", xlab="Dados simulados", ylab="Frequencia", main="Histograma")
box()
grid()
h=hist(exemplo, breaks=10, col="red", xlab="Dados simulados", ylab="Frequencia", 
       main="Histograma com a curva normal")
xfit=seq(min(exemplo),max(exemplo),length=40)
yfit=dnorm(xfit,mean=mean(exemplo),sd=sd(exemplo))
yfit=yfit*diff(h$mids[1:2])*length(exemplo)
lines(xfit, yfit, col="blue", lwd=2)
box()
grid()

Figura 1: Gráfico de histograma para dados simulados.

A ideia é mostrar que o histograma assemelha-se ao gráfico da densidade normal, a densidade dos dados.


Definição 2: Histograma, informalmente.


Matematicamente o histograma é uma função \(m_i\) que conta o número de observações que pertencem a vários intervalos disjuntos, entanto que o gráfico do histograma ou simplesmente histograma é uma mera representação desta função. Assim, se \(n\) representa o total de observações e \(k\) o número de intervalos disjuntos, o histograma satisfaz que \[ n=\sum_{i=1}^k m_i\cdot \]


O histograma é um gráfico composto por retângulos justapostos em que a base de cada um deles corresponde ao intervalo de classe e a sua altura à respectiva frequência. A construção de histogramas tem caráter preliminar em qualquer estudo e é um importante indicador da distribuição de dados. Pode indicar se uma distribuição aproxima-se de uma densidade normal como pode indicar mistura de densidades, quando os dados apresentam várias modas.

Os histogramas podem ser um mau método para determinar a forma de uma distribuição porque são fortemente influenciados pelo número de intervalos utilizados. Por exemplo, decidimos gerar 50 amostras da densidade \(\chi^2(6)\) e mostramos os gráficos de histogramas correspondentes com 14 e 26 intervalos na Figura 2.

set.seed(5678)
z=rchisq(50, df=6)
par(mar=c(5,4,2,1),mfrow=c(1,2))
hist(z, breaks=14, col="blue", main=expression(paste("Histograma ",chi^2,"(6)")), 
     ylab="Frequencia", xlab="14 intervalos")
box()
grid()
hist(z, breaks=26, col="blue", main=expression(paste("Histograma ",chi^2,"(6)")), 
     ylab="Frequencia", xlab="26 intervalos")
box()
grid()

Figura 2: Histogramas da distribuição \(\chi^2\) com 6 graus de liberdade. Número de intervalos 14 e 26, respectivamente.

Vejamos agora uma definição mais clara do histograma.


Definição 3: Histograma.


Sejam \(I_1,\cdots,I_k\) intervalos disjuntos do suporte da função de probabilidade ou de densidade da variável aleatória \(X\). O histograma é definido por \[ \widehat{f}_n(x) = \dfrac{m_i}{n|I_i|}, \quad \forall x\in I_i, \; i=1,\cdots,n, \] onde \(|I_i|\) representa o comprimento do intervalo \(i\), \(m_i\) e \(n\) como na Definição 1.


Foi provado por Robertson (1967) que, dados os intervalos \(I_1,I_2,\cdots,I_k\), o histograma \(\widehat{f}_n\) é um estimador de máxima verossimilhança dentre os estimadores expressados como funções simples e semi-contı́nuas superiormente, isto se o fecho de cada intervalos contiver duas ou mais observações. Os gráficos apresentados nas figuras acima são histogramas também segundo a proposta de Robertson (1967).

Pode-se observar que este estimador tem duas limitações importantes: a dependência do comprimento do intervalo e o fato de o histograma não constituir uma função contı́nua. A primeira destas limitações foi amplamente estudada por Wegman (1975). Ele provou que os pontos extremos de cada intervalo \(I_k\) devem ser coincidentes com observações e que, se o número mínimo de observações em cada intervalo aumenta, conforme aumenta o tamanho da amostra, o estimador \(\widehat{f}_n\) é consistente.

A segunda limitação importante do histograma, isto é, o fato de ele não constituir uma função contı́nua, incentivou diversos estudos na procura de estimadores contı́nuos da função de densidade


1.1 Cálculo automático do número de intervalos


Uma questão importante é determinar de maneira automatizada o número de intervalos disjuntos que serão utilizados para a construção do gráfico. Uma primeira forma de escolher o número de intervalos foi dada por Sturges (1926) e que constitui a forma padrão no R. Conhecida como fórmula de Sturges é dada por \[ k= [\log_2(n) + 1], \] isto significa que o número de intervalos é a parte inteira do logaritmo base 2 do número de observações mais 1.

Outras expressões comumente utilizadas são a fórmula de Scott (Scott, 1979) \(h = 3.5s/\sqrt{n}\), onde \(s\) é o desvio padrão amostral e a fórmula de Freedman-Diacconi (Freedman and Diaconis, 1981) \(h = 2\mbox{IQR}(x)/ n\), onde \(\mbox{IRQ}\) é a diferença entre o terceiro e o primeiro quantil.

Uma outra forma de definir o histograma é utilizando o estimador empírico da função de distribuição.


Definição 4: Histograma, definição moderna.


Na definição de função de densidade empírica, o parâmetro \(h\) é chamado de largura de banda. Podemos escrever o histograma \(\widehat{f}_n\), como, \[ \widehat{f}_n(x)=\dfrac{1}{2nh}\sum_{i=1}^n \pmb{1}_{(x-h;x+h)}(X_i)\cdot \]


Podemos excrever a função de densidade como \[ f(x)=\lim_{h\to 0} \dfrac{1}{h}\big(F(x+h)−F(x−h)\big), \] mas não se pode definir daqui o histograma porque esse limite é zero ou infinito e, assim, em algum momento é preciso parar, em outras palavras, não se pode chegar muito perto de zero.


Teorema 1.


Seja \(f\) a função de densidade associda à função de distribuição \(F\). Então, com probabilidade 1, \[ \widehat{f}_n(x)\sim \mbox{Binomial}(n,p), \] com \(p=F(x+h)-F(x-h)\). Assim, o comportamento assintótico do histograma pode ser derivado da distribuição binomial como \[ \mbox{E}\big(\widehat{f}_n(x)\big)=\dfrac{F(x+h)-F(x-h)}{2h} \] e \[ \mbox{Var}\big( \widehat{f}_n(x)\big)=\dfrac{p(1-p)}{2nh^2}\cdot \]


Demonstração. Exercício,


Deste teorema segue que \(\widehat{f}_n(x)\) é um estimador consistente pontual de \(f(x)\) quando \(h\to 0\) e \(nh\to\infty\). A seguir, o processo de limite é entendido como \(h=h_n\), ou seja, entende-se disto que \(h\) não será mais uma quantidade fica e sim uma quantida que depende do tamanho da amostra.

A suposição \(h=h_n\) dese satisfazer as condições do Teorema 1, de maneira que, \(h_n\to 0\) e \(nh_n\to\infty\). Estas condições podem ser interpretadas como se fosse necessário \(h_n\) ir a zero, mas não muito rápido. Isso é exatamente o que temos especulado, exceto que agora temos a taxa exata de convergência, que pode ser escrita como \(h_n=o(n)\).


Exemplo 1: Dados simulados.


Utilaremos 50 dados simulados da distribuíão \(N(0,1)\), com isso mostraremos o histograma destes 50 dados utilizando duas formas diferentes de encontrarmos uma expressão para \(h_n\), a chamada largura de banda.

set.seed(1340)
x = rnorm(50)
par(mfrow = c(1,2), mar=c(4,3,1,1))
hist(x, main = NULL, freq = FALSE, breaks = "Sturges")
box(); grid()
hist(x, main = NULL, freq = FALSE, breaks = "Scott")
box(); grid()


Neste exemplo utilizamos duas formas de escolher a largura de banda \(h_n\) dentre três diferentes possibilidades programadas na função hist. Por padrão escolhe-se breaks = “Sturges”, proposto por Sturges (1929), o qual sugere que \[ h_n=\dfrac{X_{(n)}-X_{(1)}}{1+3.322 \ln(n)}\cdot \]

A segunda situação indica que se os dados provêm da distribuição Normal temos que \(h_n=3.49sn^{−1/3}\) sendo \(s\) o desvio padrão estimado. Esta proposta deve-se à Scott (1979).

Embora o histograma é um estimador consistente quando \(h_n\to 0\) e \(nh_n\to\infty\), verifica-se que se pode fazer melhor. A melhoria também é motivada por uma preocupação prática: o histograma não é uma função suave, uma propriedade que se pode esperar que qualquer função de densidade real tenha.


2. Estimação de densidade


O estimador kernel da função de densidade pode ser realizada em qualquer número de dimensões, embora na prática a maldição da dimensionalidade faça com que seu desempenho seja degradado em dimensões altas.


Definição 5: Estimador kernel.


O estimador kernel da função de densidade é dado por \[ \widehat{f}_n(x)=\dfrac{1}{nh_n}\sum_{i=1}^n K\left(\dfrac{x-X_i}{h_n}\right), \] onde \(K(\cdot)\) é uma função conhecida como kernel.


É tipicamente assumido que \(K\) seja não-negativa, simétrica em torno de zero e satisfaz \(\int K(u)\mbox{d}u=1\). Claro que o histograma é um caso especial do estimador do kernel se \(K\) for escolhido como a função de densidade \(\mbox{Uniforme}(−1,1)\). O último não é uma função suave e é por isso que o histograma não é suave; mas escolhendo \(K\) como uma função suave, tem-se um estimador de \(f\) que seja suave.

Por exemplo, escolhendo a função de densidade \(N(0,1)\), temos por resultado o conhecido como kernel Gaussiano e assim também se utilizamos a densidade de densidade \(\mbox{Beta}\) simétrica, dada por \[ K(u)=\dfrac{\Gamma(\nu+3/2)}{\Gamma(1/2)\Gamma(\nu+1)}(1-u^2)^\nu, \quad -1<u<-1\cdot \] e \(K(u)=0\) caso contrário.

Os casos especiais \(\nu=0,1,2,3\) correspondem às funções kernel uniforme, Epanechnikov, biweight e triweight, respectivamente.

par(mfrow = c(1,2), mar=c(4,3,1,1))
kernels = eval(formals(density.default)$kernel)
plot (density(0, bw = 1), xlab = "", main = "Diferentes kernel em R"); grid()
for(i in 2:length(kernels)) lines(density(0, bw = 1, kernel =  kernels[i]), col = i)
legend(-1.0,0.2, legend = kernels, col = seq(kernels), lty = 1, cex = .8, y.intersp = 1)
h.f = sapply(kernels, function(k) density(kernel = k, give.Rkern = TRUE))
h.f = (h.f["gaussian"] / h.f)^ .2
bw = bw.SJ(x) ## escolha automática
plot(density(x, bw = bw), main = "Larguras de banda equivalentes"); grid()
for(i in 2:length(kernels)) lines(density(x, bw = bw, adjust = h.f[i], kernel = kernels[i]), col = i)
legend(55, 0.035, legend = kernels, col = seq(kernels), lty = 1)



2.1 Escolha automática da largura de banda


Um problema prático importante na estimação de densidades via kernel é como escolher a largura de banda \(h_n\) de forma automática. Note que, dadas condições como \(h_n\to 0\) e \(nh_n\to\infty\), ainda existem muitas opções para \(h_n\). Então, de certo modo, a ordem de convergência ou divergência não resolve o problema.

Uma solução para esse problema é conhecida como compensação de viés-variância. Antes de entrarmos nos detalhes, vamos primeiro apressentar um resultado em relação ao viés assintótico do estimador kernel. Aqui, o viés é definido como \[ \mbox{viés}\big( \widehat{f}_n(x)\big)=\mbox{E}\big( \widehat{f}_n(x)\big)-f(x), \] para todo \(x\in\mathbb{R}\).


Teorema 2.


Seja \(f\) uma função de densidade contínua e limitada. Então, o viés do estimador kernel de desnidade converge a zero quando \(h_n\to 0\), para todo \(x\),


Demonstração. \[ \begin{array}{rcl} \mbox{E}\big(\widehat{f}_n(x) \big) & = & \displaystyle \dfrac{1}{n}\sum_{i=1}^n \dfrac{1}{h_n} \int K\left(\dfrac{x-u}{h_n}\right)f(u)\mbox{d}u \\ & = & \displaystyle \int K(u)f(x-h_nu)\mbox{d}u \, = \, f(x) + \int K(u)\big( f(x-h_nu)-f(x)\big) \mbox{d}u\cdot \end{array} \] Utilizando então o Teorema da Convergência Dominada completa-se a demonstração.



Teorema 3.


Suponhamos que \(f\) seja contínua três vezes diferenciável, com terceira derivada limitada na vizinhança de \(x\) e \(K\), função kernel, satisfazendo que \[ \int K^2(u)\mbox{d}u< \infty \quad \mbox{e} \quad \int |u^3| K(u)\mbox{d}u < \infty \cdot \] Se \(h_n\to 0\), quando \(n\to\infty\), temos que \[ \mbox{viés}\big( \widehat{f}_n(x)\big)=\dfrac{h_n^2}{2}f''(x)\int u^2 K(u)\mbox{d}u + o(h_n^2)\cdot \] Se, além disso, \(bh_n\to\infty\) quando \(n\to\infty\), temos \[ \mbox{Var}\big(\widehat{f}_n(x)\big)=\dfrac{f(x)}{nh_n}\int K^2(u)\mbox{d}u+o(1/nh_n)\cdot \]


Demonstração. A demonstração é baseada na expansão de Taylor, \[ f(x-h_n u)=f(x)-h_n u f' (x)+\dfrac{h_n^2 u^2}{2} f''(x)-\dfrac{h_n^3 u^3}{6}f''' (\epsilon), \] sendo que \(\epsilon\) fica entre \(x−h_n u\) e \(x\). Os detalhes são deixados como um exercício.


Uma medida de precisão do estimador é o erro quadrático médio (EQM), dado por \[ \mbox{EQM}\big( \widehat{f}_n(x)\big)=\mbox{E}\big(\widehat{f}_n(x)-f(x)\big)^2\cdot \] É fácil mostrar que o EQM combina o viés e a variância de tal maneira que \[ \mbox{EQM}\big( \widehat{f}_n(x)\big)=\mbox{viés}\big(\widehat{f}_n(x)\big)^2+\mbox{Var}\big(\widehat{f}_n(x)\big)\cdot \]

Vemos que, sob as condições \(h_n\to 0\) e \(nh_n\to\infty\) e, se ignorarmos os termos de baixa ordem, temos \[ \mbox{EQM}\big( \widehat{f}_n(x)\big) \approx \dfrac{h_n^4}{4}\big( f''(x)\big)^2\tau^4+\dfrac{f(x)}{nh_n}\gamma^2, \] onde \(\tau^2=\int u^2 K(u)\mbox{d}u\) e \(\gamma^2=∫ K^2(u)\mbox{d}u\).

O termo à direita da expressão acima é minimizada quando \[ h_n=\left( \dfrac{\gamma^2 f(x)}{\tau^4 \big( f''(x)\big)^2}\right)^{1/5}n^{-1/5}\cdot \] Note ainda que a expressão acima não é a solução ideal, isso porque \(f\) é desconhecida na prática. No entanto, dá-nos pelo menos alguma ideia sobre a taxa ideal de convergência a zero, sendo esta \(h_n=O(n^{-1/5})\).

Quando \(f\) é desconhecida, uma abordagem natural seria substituí-lo por um estimador e, assim, obter uma largura de banda ideal estimada. Uma complicação é que a largura de banda ideal depende de \(x\) mas, idealmente, gostaríamos de usar uma largura de banda que funcionasse para diferentes \(x\) dentro de um certo intervalo, se não todos os \(x\).

Para obter uma largura de banda ideal que não depende de \(x\), integramos os dois lados da expressão do EQM em relação a \(x\). Isto nos leva a \[ \int \mbox{EQM}\big(\widehat{f}_n(x)\big)\approx \dfrac{\tau^4 h_n^4}{4}\int \big(f''(x)\big)^2\mbox{d}x+\dfrac{\gamma^2}{nh_n}\int f(x)\mbox{d}x = \dfrac{\tau^4 \theta^2 h_n^4}{4}+\dfrac{\gamma^2}{nh_n}, \] com \(\theta^2=\displaystyle \int \big( f''(x)\big)^2 \mbox{d}x\).

Pelo mesmo argumento, o lado direito acima é minimizado quando \[ h_n=\left( \dfrac{\gamma^2}{\tau^4\theta^2}\right)^{1/5}n^{-1/5}\cdot \]

Desta vez, o \(h_n\) ideal não depende de \(x\). Além disso, a integral do EQM ou o IEQM mínimo é dado por \[ \mbox{IEQM}=\int \mbox{EQM}\big(\widehat{f}_n(x)\big)\mbox{d}x=\dfrac{5}{4}\big(\tau\gamma^2\big)^{4/5}\theta^{2/5}n^{-4/5}\cdot \]

Uma implicação do resultado acima é a seguinte. Note que o IEQM depende do kernel \(K\) através de \(c_K=\big(\tau\gamma^2\big)^{4/5}\). Mostrou-se que para os kernels comumente usados, tais como aqueles listados, o desempenho dos estimadores de kernel correspondentes é quase o mesmo em termos dos valores de \(c_K\).

Voltando ao problema sobre a estimação da largura de banda ideal, vemos que tudo o que precisamos é encontrar um estimador consistente de \(\theta^2\). Se \(f\) é a função de densidade da distribuição normal com desvio padrão \(\sigma\), então pode ser mostrado que \(\theta^2=3/8\sqrt{\pi}\sigma^5\). Naturalmente, se alguém souber que \(f\) é normal, então a estimação da densidade não-paramétrica não seria necessária, porque um método paramétrico provavelmente seria melhor. Em geral, pode-se expandir \(f\) em torno da densidade gaussiana usando a expansão de Edgeworth.

Utilizando a abordagem acima, Hjort and Jones (1996) obteveram o seguinte estimador ótimo para a largura de banda \[ \widehat{h}_n=\widehat{h}_0\left( 1+\dfrac{35}{48}\widehat{\gamma}_4+\dfrac{35}{32}\widehat{\gamma}^2_3+\dfrac{385}{1024}\widehat{\gamma}^2_4\right)^{-1/2}, \] onde \(\widehat{h}_0\) é a estimativa ideal da largura de banda assumindo que \(f\) é normal, isto é, com \(\theta^2\) substituído por \(3/8\sqrt{\pi}\sigma^5\) ou mais explicitamente \[ \widehat{h}_0=1.06\left( \dfrac{\widehat{\sigma}}{n^{1/5}}\right)\cdot \]

Chamamos \(\widehat{h}_0\), a largura de banda da linha de base e \(\widehat{\sigma}^2\), a variância amostral, dada por \[ \widehat{\sigma}^2 = \dfrac{1}{n-1}\sum_{i=1}^n (x_i-\overline{x})^2\cdot \]

Além disso, \(\widehat{\gamma}_3\) e \(\widehat{\gamma}_4\) são os estimadores dos coefcientes de assimetria e curtose amostrais, respectivamente, dados por \[ \widehat{\gamma}_3=\dfrac{1}{(n-1)\widehat{\sigma}^3}\sum_{i=1}^n (x_i-\overline{x})^3 \] e \[ \widehat{\gamma}_4=\dfrac{1}{(n-1)\widehat{\sigma}^4}\sum_{i=1}^n (x_i-\overline{x})^4-3\cdot \]

Os procedimentos de estimação de densidade mostrados até o momento foram implementados na função density.


Exemplo 2: Estimador kernel de densidade em dados simulados.


Utilaremos os dados simulados da distribuíão \(N(0,1)\) no exemplo anterior. Mostraremos o histograma e o estimador kernel da função de densidade.

par(mfrow = c(1,2), mar=c(4,3,1,1))
hist(x, main = NULL, freq = FALSE, breaks = "Sturges")
box(); grid()
densidade1 = density(x, bw = "nrd0")
lines(densidade1, col = "red")
par(mar=c(4,2,1,1))
hist(x, main = NULL, freq = FALSE, breaks = "Scott")
box(); grid()
densidade2 = density(x, bw = "bcv")
lines(densidade2, col = "red")


Houveram outras abordagens para a seleção da largura de banda ótima, incluindo o método de validação cruzada.


2.2 Validação cruzada


No contexto da estimativa da densidade do kermel e da seleção da largura de banda, a largura de banda não deve ser muito pequena porque a variância seria muito grande e não deve ser muito grande, o viés seria muito grande. Outro método padrão para selecionar a largura de banda é a técnica de validação cruzada, descrita em Chiu (1991).

Aqui, gostaríamos de minimizar \[ \mbox{E}\left(\int\Big(\widehat{f}_n(u)-f(u) \Big)^2\mbox{d}u \right)\cdot \] A integral pode ser escrita \[ \int \widehat{f}_n^2(u) \mbox{d}u -2\int \widehat{f}_n^2(u)f(u)\mbox{d}u+\int f²(u)\mbox{d}u\cdot \]

Como a terceira componente é constante, temos que minimizar o valor esperado da soma das duas primeiras.

A ideia é aproximá-lo como \[ J(h)=\int \widehat{f}_n^2(u) \mbox{d}u -\dfrac{2}{n}\sum_{i=1}^n \widehat{f}_{n_{(-i)}}^2(x_i) \] que pode ser facilmente calculado onde \(\widehat{f}_{n_{(-i)}}\) é o estimador de densidade obtido após a remoção da \(i\)-ésima observação. Referimo-nos a \(J(h)\) como o escore de validação cruzada ou risco estimado.


Exemplo 3: Estimador kernel por validação cruzada.


O objetivo deste exemplo é comparar os valores das larguras de banda \(\widehat{h}_n\) obtida pelo método da integral do EQM e àquela obtida por validação cruzada.

names(densidade1)
## [1] "x"          "y"          "bw"         "n"          "old.coords"
## [6] "call"       "data.name"  "has.na"
densidade1$bw
## [1] 0.3635834
densidade2$bw
## [1] 0.5098764


Otimizando a função \(J(h)\) encontramos o ponto mínimo e o valor mínimo como a seguir.

J=function(h){
  fhat=Vectorize(function(y) density(x,from=y,to=y,n=1,bw=h)$y)
  fhati=Vectorize(function(i) density(x[-i],from=x[i],to=x[i],n=1,bw=h)$y)
  F=fhati(1:length(x))
  return(integrate(function(y) fhat(y)^2,-Inf,Inf)$value-2*mean(F))
}
vx=seq(0.1,0.9,by=.01)
vy=Vectorize(J)(vx)
df=data.frame(vx,vy)
library(ggplot2)
qplot(vx,vy,geom="line",data=df)


optimize(J,interval=c(.1,0.9))
## $minimum
## [1] 0.606255
## 
## $objective
## [1] -0.2637266

Observe que, de fato, está próximo da largura de banda ideal de Silverman.



3. Exemplo de aplicação


Ilustraremos a suavização usando um pequeno conjunto de dados sobre ângulos formados entre dois planos peptídicos subsequentes em estruturas proteicas 3D. Este conjunto de dados é selecionado porque as distribuições angulares são multimodais e ligeiramente não padronizadas e essas propriedades são adequadas para ilustrar considerações fundamentais sobre a estimativa de densidade suave na prática.

Figura A: Estrutura 3D das proteínas é amplamente dada pelos ângulos \(\phi\) e \(\psi\), dos planos peptídicos. Imagem retirada do livro “Computational Statistics with R”, disponível no endereço: https://cswr.nrhstat.org/intro.

Consideraremos um pequeno conjunto de dados phipsi, de ângulos determinados experimentalmente a partir de uma única proteína, a proteína humana 1HMP, que é composta por duas cadeias, denotadas A e B. Mostra-se na Figura A a estrutura 3D da proteína.

Uma proteína é uma molécula grande que consiste em uma espinha dorsal com átomos de carbono e nitrogênio dispostos sequencialmente: Figura B: Espinha dorsal de uma proteina.

Um átomo de hidrogênio se liga a cada nitrogênio (N) e um átomo de oxigênio se liga a cada carbono sem o subscrito \(\alpha\) (C), veja a figura acima, e esses quatro átomos formam juntos o que é conhecido como ligação peptídica entre dois átomos de carbono alfa (\(\mbox{C}_\alpha\)).

Cada átomo de \(\mbox{C}_\alpha\) liga-se a um átomo de hidrogênio e a uma cadeia lateral de aminoácido. Existem 20 aminoácidos naturais em proteínas geneticamente codificadas, cada uma com um código de três letras, como Gly para Glicina, Pro para Prolina, etc.. A proteína normalmente formará uma estrutura 3D complicada determinada pelos aminoácidos, que por sua vez determinam ângulos \(\phi\) e \(\psi\) entre os planos peptídicos como mostrado na Figura A.

phipsi = read.csv("http://leg.ufpr.br/~lucambio/Nonparam/phipsi.csv", sep = ",", header = T)
head(phipsi)
##   chain  AA pos        phi        psi
## 1     A Pro   5  -92.92684   12.94131
## 2     A Gly   6   65.79681 -162.22971
## 3     A Val   7  -81.13208  121.41302
## 4     A Val   8  -85.52381  137.17372
## 5     A Ile   9 -124.98875   85.24247
## 6     A Ser  10  -43.11722  147.11371


Alguma estatísticas descritivas da variável \(\psi\).

summary(phipsi$psi)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -178.517  -37.730    9.247   40.277  133.481  178.871


Transformando os dados da escala em graus para a escala em radianos.

deg2rad = function(deg) {(deg * pi) / (180)}
summary(deg2rad(phipsi$psi))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -3.1157 -0.6585  0.1614  0.7030  2.3297  3.1219


Podemos usar funções de base R, como hist e density, para visualizar as distribuições marginais dos dois ângulos.

library(dplyr)
par(mar=c(4,4,1,1), mfrow=c(1,2))
hist(deg2rad(phipsi$psi), prob = TRUE, ylab="Densidade", xlab = expression(phi),
     main=expression(paste("Histograma de ",phi)))
rug(deg2rad(phipsi$psi))
grid()
density(deg2rad(phipsi$psi)) %>% lines(col = "red", lwd = 2)
box()
hist(phipsi$psi, prob = TRUE, ylim=c(0,0.008), ylab="Densidade", xlab = expression(psi),
     main=expression(paste("Histograma de ",psi)))
rug(phipsi$psi)
density(phipsi$psi) %>% lines(col = "red", lwd = 2)
box()
grid()


A curva suave de densidade vermelha mostrada na figura acima pode ser considerada uma versão suave de um histograma. É surpreendentemente difícil encontrar procedimentos de suavização automática que funcionem uniformemente bem – é até bastante difícil selecionar automaticamente o número e as posições das quebras usadas para histogramas. Este é um dos pontos importantes abordados na literatura: como implementar boas escolhas padrão de vários parâmetros de ajuste que são exigidos por qualquer procedimento de suavização.

Também é possível usar o pocote ggplot2 para obter resultados semelhantes.

library(ggplot2)
grafico1 = ggplot(phipsi, aes(x = phi)) + 
  geom_histogram(aes(y = after_stat(density)), bins = 13) + 
  geom_density(col = "red", linewidth = 1) + 
  geom_rug()
grafico2= ggplot(phipsi, aes(x = psi)) + 
  geom_histogram(aes(y = after_stat(density)), bins = 13) + 
  geom_density(col = "red", linewidth = 1) + 
  geom_rug()
library(gridExtra)
grid.arrange(grafico1,grafico2, widths=c(1,2))


O intervalo dos dados do ângulo é conhecido como \((-\pi,\pi]\), do qual nem o histograma nem o suavizador de densidade aproveitam. Observe também que para o ângulo \(\psi\), parece que os padrões resultam em suavização excessiva da estimativa de densidade. Ou seja, a densidade é mais suavizada do que os dados e o histograma parecem suportar.

Para obter resultados diferentes – e talvez melhores – podemos tentar alterar alguns dos padrões do histograma e das funções de densidade. Os dois padrões mais importantes a serem considerados são a largura de banda e o kernel. O kernel controla como os pontos de dados vizinhos são ponderados relativamente uns aos outros e a largura de banda controla o tamanho das vizinhanças.

Uma largura de banda pode ser especificada manualmente como um valor numérico específico, mas para um procedimento totalmente automático, ela é selecionada por um algoritmo de seleção de largura de banda. O padrão density é um algoritmo bastante simplista conhecido como regra prática de Silverman.

hist(deg2rad(phipsi$psi), breaks = seq(-pi-0.1, pi+0.1, length.out = 15), prob = TRUE, 
     ylab="Densidade", xlab = expression(phi), main=expression(paste("Histograma de ",phi," em graus")))
rug(phipsi$psi)
density(deg2rad(phipsi$psi), adjust = 1, cut = 0) %>% lines(col = "red", lwd = 2)
density(deg2rad(phipsi$psi), adjust = 0.5, cut = 0) %>% lines(col = "blue", lwd = 2)
density(deg2rad(phipsi$psi), adjust = 2, cut = 0) %>% lines(col = "purple", lwd = 2)
box()
grid()

hist(deg2rad(phipsi$phi), breaks = seq(-pi-0.1, pi+0.1, length.out = 15), prob = TRUE, ylab="Densidade", 
     xlab = expression(psi), main=expression(paste("Histograma de ",psi," em graus")))
rug(phipsi$phi)
# Default kernel is "gaussian"
density(deg2rad(phipsi$phi), bw = "SJ", cut = 0) %>% lines(col = "red", lwd = 2)
density(deg2rad(phipsi$phi), kernel = "epanechnikov", bw = "SJ", cut = 0) %>% 
  lines(col = "blue", lwd = 2)
density(deg2rad(phipsi$phi), kernel = "rectangular", bw = "SJ", cut = 0) %>% 
  lines(col = "purple", lwd = 2)
box()
grid()

Estes gráficos mostram os histogramas e várias estimativas de densidade para os âgulos \(\phi\) e \(\psi\). As cores indicam diferentes opções de ajustes de largura de banda usando a seleção de largura de banda padrão e diferentes opções de kernels usando a seleção de largura de banda Sheather-Jones (SJ).

Os gráficos mostram exemplos de diversas estimativas de densidade diferentes que podem ser obtidas alterando os padrões de densidade. As quebras do histograma também foram escolhidas manualmente para garantir que correspondam ao intervalo dos dados. Observe, em particular, que a seleção da largura de banda Sheather-Jones (SJ) parece funcionar melhor do que a largura de banda padrão deste exemplo. Este é geralmente o caso de distribuições multimodais, onde a largura de banda padrão tende a suavizar demais. Observe também que a escolha da largura de banda tem muito mais consequências do que a escolha do kernel, esta última afetando principalmente o quão instável é a estimativa de densidade localmente.

Deve-se notar que os padrões surgem como uma combinação de escolhas historicamente sensatas e compatibilidade retroativa. Deve-se pensar na escolha de um padrão bom e robusto, mas uma vez escolhido um padrão, ele não deve ser alterado ao acaso, pois isso pode quebrar o código existente. É por isso que nem todos os padrões usados em R são, pelos padrões atuais, as escolhas mais conhecidas. Você vê esse argumento na documentação de density em relação ao padrão para seleção de largura de banda, onde Sheather-Jones (SJ) é sugerido como um padrão melhor que o atual, mas por razões de compatibilidade, a regra prática de Silverman é o padrão e provavelmente permanecerá sendo então.


4. Avaliação do número de modas


Dada uma amostra de dados de uma variável aleatória, determinar o número de modas na densidade é uma questão relevante para apoiar futuras decisões durante a abordagem de modelagem. É claro que distribuições unimodais, como a densidade gaussiana, podem não ser adequadas para caracterizar o comportamento de mecanismos geradores de dados mais complexos em ciências aplicadas.

Alguns exemplos que requerem distribuições mais complexas para refletir o número real de modas podem ser encontrados em muitos campos aplicados, como a astronomia, por exemplo, no estudo de padrões unimodais ou multimodais dos períodos de rotação de estrelas para diferentes temperaturas; administração de empresas, por exemplo, ao analisar o capital investido em campanhas de crowdfunding; ciência florestal, por exemplo, na análise do número de modas na distribuição de medidas de retroespalhamento, para áreas de floresta densa e sem vegetação, dependendo da porcentagem de pixels no solo; genética, por exemplo, para identificar quais CpGs, regiões de DNA onde um nucleotídeo de citosina é seguido por um nucleotídeo guanina, apresentam distribuições multimodais ou psicologia, onde, por exemplo, o estudo do número das modas é crucial para detectar a presença de fenômenos cognitivos de processo único ou duplo; entre outros.

Em princípio, suavizadores de densidade não paramétricos, como o estimador de densidade kernel, introduzido por Rosenblatt (1956) e Parzen (1962), podem superar o problema de restringir a estimativa da função de densidade a uma família paramétrica previamente especificada. No entanto, surgem duas questões importantes ao realizar a estimativa da densidade por meio de métodos de suavização de densidade de kernel ou qualquer outro. A primeira questão é que os profissionais podem sentir-se mais confortáveis ao interpretar e lidar com modelos paramétricos, uma vez que em muitos casos as estimativas dos parâmetros podem ser interpretadas em termos da distribuição dos dados, uma vez que controlam algumas características específicas.

A segunda questão é que, mesmo estando satisfeito com o resultado não paramétrico do estimador de densidade kernel, uma vez que fornece uma versão estimada da distribuição, pode haver algumas dúvidas sobre as características destacadas por este estimador: são genuínas a partir da distribuição ou são apenas devido à variabilidade da amostragem?

As preocupações anteriores podem ser parcialmente resolvidas ou respondidas pela identificação das modas significativas no estimador de densidade kernel. Assim, como passo anterior antes de ajustar um modelo paramétrico, deve-se verificar quantos grupos distinguíveis existem na distribuição dos dados, sendo estes grupos identificados pelas modas da densidade.

Isto pode ser feito por métodos exploratórios ou por procedimentos de teste. Em ambos os casos, deve-se determinar também quanto do padrão observado no estimador de densidade é real e quanto se deve a artefatos de amostragem. Além disso, uma aproximação paramétrica muito flexível, mas simples, com vários grupos/modos pode ser realizada ajustando misturas de normais.

Ferramentas exploratórias ou descritivas, embora possam forner uma análise completa da estimativa de densidade, requerem uma decisão sobre o número de modas a serem considerados após examinar uma saída gráfica. Portanto, as conclusões não podem ser obtidas diretamente através da aplicação de um procedimento automático que indique quantas das modas observados nas representações gráficas são realmente significativas.

Esta questão pode ser respondida por um teste de hipótese: \(H_0 : j = k\) vs. \(H_A : j > k\), denotando por \(j\) o número real de modas na densidade e sendo \(k\) um número inteiro positivo, então \(k = 1\) é um teste de unimodalidade. Este problema de teste foi resolvido projetando estatísticas de teste baseadas na largura de banda crítica e/ou no excesso de massa. Esses procedimentos serão brevemente descritos, juntamente com métodos exploratórios.


4.1 Ferramentas exploratórias para avaliar a multimodalidade


Como o número de modas em \(\widehat{f}_{h_n}\) é uma função decrescente monótona de \(h_n\), quando o kernel gaussiano é usado, uma solução exploratória simples, para determinar o número de modas, representa esta estimativa de densidade para diferentes valores de \(h_n\). Na verdade, esta é a ideia em algumas ferramentas gráficas, como a árvore de modas e a floresta de modos, onde um exemplo de ambas as representações é fornecido na Figura 1 (painéis (b) e (c)).

Na árvore de modos, Minnotte e Scott (1993) criaram um diagrama de árvore (semelhante ao dendograma) representando, com linhas verticais contínuas, as localizações dos modos (eixo primário) de ˆfh para diferentes parâmetros de largura de banda h (eixo secundário). Além disso, representa, com linhas tracejadas horizontais, como cada modo se divide em mais modos à medida que a largura de banda diminui (de cima para baixo), mostrando a relação entre os novos modos e os modos originais dos quais eles se dividiram.

Como apontado por Minnotte et al. (1998), o problema da árvore modal é a forte dependência da amostra disponível. Essa é a razão pela qual a floresta modal é construída calculando a posição dos modos estimados a partir de diferentes árvores modais obtidas a partir da amostragem com substituição da amostra original. Para facilitar a visualização deste método exploratório, a janela gráfica é dividida em diferentes pixels (previamente escolhidos) de localização-largura de banda (eixo horizontal-vertical). Então, esta ferramenta representa o número de vezes que um modo estimado cai em cada pixel (largura de banda de localização), sombreando-o proporcionalmente às contagens (contagens grandes correspondendo a pixels mais escuros e contagens baixas a pixels mais claros). Então, na floresta de modos, os modos são identificados por regiões cinza escuro.

Esta seção fornece um breve histórico sobre diferentes ferramentas exploratórias e de testes incluídas no multimodo. Um elemento chave nos fundamentos das diferentes propostas é o estimador de densidade kernel. Dada uma amostra aleatória \(X_1,\cdots, X_n\) de uma variável aleatória \(X\) com densidade desconhecida \(f\), o estimador de densidade kernel para um x ∈ R fixo é definido como

library("multimode")
data("stamps", package = "multimode")
hist(stamps, prob = TRUE, ylab="Densidade", xlab = "Espessura", main="Histograma da espessura dos selos")
rug(stamps)
density(stamps) %>% lines(col = "red", lwd = 2)
box()
grid()

bwRT = bw.nrd0(stamps)
bwRT
## [1] 0.003909682
bwPI = bw.SJ(stamps)
bwPI
## [1] 0.001205108
nmodes(data = stamps, bw = bwRT)
## [1] 2
nmodes(data = stamps, bw = bwPI, full.result = TRUE)
## 
## n= 485. Number of modes: 9
## Bandwidth: 0.001205108
## Support where the number of modes is computed
##      -Inf      Inf


O conjunto de dados de stamps (selos), analisados em Izenman and Sommer (1988), consiste em medições de espessura (em milímetros) de 485 selos de tecido branco usados, sem marca d’água, da edição de 1872 do selo Hidalgo do México. Todos tinham sobreimpressão com o ano 1872 ou 1873 ou 1874 e alguns deles tinham marca d’água (Papel Sellado ou LA+-F ), estando esta informação também incluída no interior dos selos.

Dado que o valor dos selos depende da sua escassez, é importante determinar o número de grupos disponíveis numa determinada emissão de selos. Para esta questão específica de selos, embora a marca de água em alguns selos (em 29 de 485) ajude a concluir que existem pelo menos dois grupos, a questão sobre o número de grupos pode ser respondida analisando o número subjacente de modas.


5. Estimação de densidade por polinômios de Legendre


A forma de estimar a função de densidade descrita nesta seção foi proposta por Tien-Wen (1979). Descrevemos aqui a teoria envolvida para entendermos a proposta, assim como sua implementação no R.

Consideramos \(X_1,\cdots,X_n\) uma amostra aleatória, ou seja, um conjunto de variáveis independentes igualmente distribuídas com densidade desconhecida \(f\). Queremos estimar \(f(x)\), \(x\in\mathbb{R}\), utilizando polinômios de Legendre. Estes polinômios foram definidos por Adrien-Marie Legendre, um matemático francês, por volta do final do século XVIII e início do século XIX.


Definição ?: Polinômios de Legendre.


Seja \(-1\leq x\leq 1\), dizemos que \(P_n(x)\) é o \(n\)-ésimo polinômio de Legendre se \[ P_n(x)=\dfrac{1}{2^n n!}\dfrac{\mbox{d}^n}{\mbox{d}x^n} (x^2-1)^n\cdot \]


Uma propriedade importante, demonstrado por Sansone (1959), é a seguinte \[ |P_n(x)|\leq 1, \qquad n=0,1,2,\cdots \] Usando o teorema binomial em \((x^2-1)^n\) e diferenciando o resultado \(n\) vezes, obtemos que \[ P_n(x)=\sum_{k=0}^N (-1)^k \dfrac{(2n-2k)!}{2^n k! (n-k)! (n-2k)!} x^{n-2k}, \] onde \(N=n/2\) se \(n\) for par e \(N=(n-1)/2\) caso \(n\) for ímpar.

A figura abaixo mostra os primeiros 6 polinômios de Legendre, para isso utilizamos a função legendre.polynomials no pacote orthopolynom. Esta mesma função, mas agora acresentando uma opção e escrevendo legendre.polynomials(n, normalized=FALSE), devolve como valor uma lista com os polinômios de Legendre de 0 a \(n\). Esta lista de caracteres deve ser transformada numa expressão e para isso utilizamos a função parse e então utilizamos eval para avaliar essa expressão como uma função.

library(orthopolynom)
# Cria uma lista com os polinômios de Legendre de grau 0 a 5
Legendre <- legendre.polynomials(5, normalized=FALSE)
Legendre
## [[1]]
## 1 
## 
## [[2]]
## x 
## 
## [[3]]
## -0.5 + 1.5*x^2 
## 
## [[4]]
## -1.5*x + 2.5*x^3 
## 
## [[5]]
## 0.375 - 3.75*x^2 + 4.375*x^4 
## 
## [[6]]
## 1.875*x - 8.75*x^3 + 7.875*x^5
# Convertendo uma cadena de caracteres, resposta da função legendre.polynomials,  
# em uma expressão
x <- seq(-1, 1, length.out = 100)
legendre <- matrix(0,nrow=length(x),ncol=6)
for (i in 1:6){
  funcao_string <- paste("function(x)",Legendre[[i]])
  funcao <- eval(parse(text = funcao_string))
  legendre[,i] <- funcao(x) 
}
# Plotar os polinômios de Legendre
par(mar=c(4,4,1,1))
plot(x, legendre[,1], type = "l", col = "red", ylim = c(-1, 1), 
     main = "Polinômios de Legendre de Grau 0 a 5", xlab = "x", ylab = "P_n(x)")
for (i in 2:6) {
  lines(x, legendre[,i], col = rainbow(6)[i])
}
legend("bottomright", legend = 0:5, col = rainbow(6), lty = 1, bty = "n")
grid()

Figura ?: Polinômios de Legendre de ordem 0 a 5.

Assumindo que \(f(x)\in L^2[-1,1]\), sendo denotado por \(L^2[-1,1]\) o conjunto de todas as funções reais de quadrado integráveis no intervalo fechado [-1,1], foi demonstrado por Hardy and Rogosinski (1950) que a função de densidade pode ser escrita como \[\begin{equation} \tag{5.1} f(x)=\sum_{n=0}^\infty a_n P_n(x)\sqrt{(2n+1)/2}, \end{equation}\] sendo os coneficientes \(a_n\) definidos como \[\begin{equation} \tag{5.3} a_n = \sqrt{(2n+1)/2} \int_{-1}^1 f(x) P_n(x)\mbox{d} x\cdot \end{equation}\]

Podemos perceber que, a expressão em (5.1) ainda não é útil em situações práticas por dois motivos: primeiro é uma construção com infinitos termos, devemos parar a soma numa quantidade finita de termos mas, em quantos? qual a qualidade dessa aproximação por finitos termos? em segundo lugar, os coeficientes \(a_n\) dependem da função de densidade dos dadps \(f\), desconhecida por nós.

a <- rep(0,6)
a[1] = sqrt(1/2)*eval(Legendre[[1]])*integrate(dnorm, lower = -1, upper = 1)$value
for (i in 2:6){
  funcao_string <- paste("function(x) (", Legendre[[i]],")*exp(-x^2/2)")
  funcao <- eval(parse(text = funcao_string))
  a[i] = sqrt((2*(i-1)+1)/2)*sqrt(1/(2*pi))*integrate(funcao, lower = -1, upper = 1)$value
}
print(a, digits = 5)
## [1]  0.4827344  0.0000000 -0.0683411  0.0000000  0.0043012  0.0000000
par(mar=c(4,4,1,1),pch=19)
curve(dnorm(x), from = -3, to = 3, col = "red", xlab = "X ~ N(0,1)", ylab = "")
rug(a)

x <- seq(-3,3,by=0.01)
#density.Legendre=function(x){
#  f=x; n=length(x)
#  for (i in 1:n){
#    for (j in 1:6){
#      f[i,j] = a[j]*Legendre[j]*sqrt((2*(j-1)+1)/2)
#    }
#  }
#  for (i in 1:n){
#    g[i] = sapply(f,)
#  }
#  f[i] = a[i]Vectorize(function(y) density(x,from=y,to=y,n=1,bw=h)$y)
#  fhati=Vectorize(function(i) density(x[-i],from=x[i],to=x[i],n=1,bw=h)$y)
#  F=fhati(1:length(x))
#  return(integrate(function(y) fhat(y)^2,-Inf,Inf)$value-2*mean(F))
#  
#}
densidade.estimada = a[1]*eval(Legendre[[1]])
#curve(densidade.Legendre(rnorm(30)), col="blue")
#grid()


Acontece que podemos escrever de uma outra maneira a função de densidade como \[\begin{equation} \tag{5.2} f(x)=\sum_{n=0}^\infty a_n \psi_n(x), \end{equation}\] sendo que \[ a_n = \int_{-1}^1 f(x)\psi_n(x)\mbox{d} x=\int_{-1}^1 f(x)\sqrt{(2n+1)/2}P_n(x)\mbox{d} x\cdot \]

Na definição em (1) utilizamos as chamadas funções normalizadas de Legendre, definidas como \[ \psi_n(x)=\left[\dfrac{1}{2}(2n+1)\right]^{1/n} P_n(x), \qquad n=0,1,2,,\cdot \] Se \(L\) Estas funções formam um sistema complto ortogonal

\[ \widehat{f}_n(x)=\sum_{k=0}^{l(n)} \widehat{a}_{kn} \psi_k(x), \] com \[ \widehat{a}_{kn}=\dfrac{1}{n}\sum_{j=1}^n \psi_k(x_j), \] sendo que \(\psi_k(\cdot)\) é a \(k\)-ésima função de Legendre padronizada e \(l(n)\) é um inteiro que depende do tamanho da amostra \(n\).

Colocarmos algumas suposições em \(f\) podemos demonstrar que \(\widehat{f}_n(x)\) é consistente no sentido do erro quadrático médio integral e erro quadrático médio. Podemos também encontrar a taxa de convergência do estimador à \(f(x)\).


6. Bibliografia


Hardy, G. H., and W. W. Rogosinski. 1950. Fourier Series. Cambridge University Press.
Robertson, T. 1967. “On Estimating a Density Which Is Measurable with Respect to a σ-Lattice.” The Annals of Mathematical Statistics, no. Vol.38 (2): 482–93.
Sansone, G. 1959. Orthogonal Functions. Interscience, New York.
Sturges, H. 1926. “The Choice of a Class-Interval.” Journal of the American Statistical Association, no. 21: 65–66.
Tien-Wen, C. 1979. “On the Probability Density Estimation by Legendre Polynomial.” Soochow Journal of Mathematics, no. 5: 163–69.
Wegman, E. J. 1975. “Maximum Likelihood Estimation of a Probability Density Function.” Sankhyā: The Indian Journal of Statistics, Series A, no. 37(2) http://www.jstor.org/stable/25049976: 211–24.