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 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çã de \(\theta\). Em um problema não-paramétrico, a distribuição é totalmente desconhecida com, talvez, algumas restriçõs 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"         "call"      "data.name"
## [7] "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.6083966
## 
## $objective
## [1] -0.2637472

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 subjacente é 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 (McQuillan, Mazeh and Aigrain 2014); administração de empresas, por exemplo, ao analisar o capital investido em campanhas de crowdfunding (Colombo, Franzoni and Rossi-Lamastra 2015); 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 (Santoro et al. 2011); 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 (Joubert et al. 2016); 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 (Freeman and Dale 2013); 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 de densidade a uma família paramétrica previamente especificada. No entanto, surgem duas questões importantes ao realizar a estimativa de 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 subjacente, pode haver algumas dúvidas sobre as características destacadas por este estimador de curva: 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, uma revisão sobre este tópico pode ser encontrada, por exemplo, em McLachlan and Peel (2000).

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.