Moda de uma amostra: método do histograma vs método kernel

November 18, 2011
By

This post was kindly contributed by Ridículas - go there to comment and to read the full post.

Moda de uma amostra aleatória de uma distribuição contínua obtida pelo método do histograma e da densidade kernel.

Uma pergunta bem comum nos fóruns de estatística, geralmente associado à alunos de graduação, é de alguma função no R para obter a moda de um conjunto de dados. Bem, como moda se refere ao valor de maior frequência na amostra, a tarefa é simples para dados discretos: basta obter a distribuição de frequência dos valores (função table()) e tomar aquele associado à maior frequência.

Para dados contínuos, teoricamente, não teríamos dois valores observados iguais em uma amostra. Em amostras reais às vezes temos, seja pelo fato da amostra ser grande e/ou pelo instrumento de medida não ser preciso para diferenciar duas observações. Por exemplo, uma balança que registra massa em quilos vai estar sempre arredondando (censurando) a parte decimal dos valores. Então na prática até que podemos ter, eventualmente, valores repetidos em amostras de variáveis contínuas. Porém, não é assim que calculamos moda nessas amostras, porque valores repetidos podem não ocorrer e quando ocorrem, raramente excede à duas ocorrência de um mesmo valor.

A obtenção da moda pelo método do histograma é simples e está bem documentada nos livros texto da disciplina de estatística básica. Consiste em agrupar os dados em classes e fazer uma interpolação linear que considera a frequência da classe modal e de suas vizinhas. A interpolação linear usa semelhança de triângulos para encontrar a moda da amostra. A fórmula é essa

mo= l_i+\displaystyle A_c \left(\frac{f_m-f_e}{2f_m-f_e-f_d}\right )

em que l_i é o limite inferior da classe modal, A_c é a amplitude da classe modal, f_m é a frequência da classe modal, f_e é a frequência da classe à esquerda da modal, f_d é a frequência da classe à direita da modal.

Outra forma de obter a classe modal é calculando o valor da amostra que corresponde ao máximo da função de densidade kernel. É um método mais moderno e computacionalmente mais trabalhoso, pois deve-se ajustar a densidade kernel, fazer aproximação da função kernel por uma função interpoladora e depois obter o máximo dessa função. Mas como somos usuários de R, já temos há um bom tempo as ferramentas necessárias para isso.

No CMR abaixo eu implementei os métodos para obtenção da moda de uma amostra pelos dois métodos: do histograma e da densidade kernel. Até a próxima ridícula.

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição discreta

x <- rpois(100, lambda=12)
plot(sort(x))
tb <- table(x); tb
tb[which.max(tb)]

#-----------------------------------------------------------------------------
# amostra aleatória de distribuição contínua

x <- rbeta(1000, 0.5,0.5)
ht <- hist(x) # histograma

#-----------------------------------------------------------------------------
# função para obter a moda pelo método do histograma, moda obtida por
# interpolação linear, baseado em semelhaça de triângulos

moda.hist <- function(ht, plotit=TRUE){
  ## ht: um objeto do uso da função hist()
  mcl <- which.max(ht$counts)
  li <- ht$breaks[mcl]
  width <- diff(ht$breaks[mcl+0:1])
  counts <- c(0,ht$counts,0)
  delta <- abs(diff(counts[1+mcl+(-1:1)]))
  moda <- li+width*delta[1]/sum(delta)
  cols <- rep(5, length(ht$counts))
  cols[mcl] <- 3
  if(plotit==TRUE){
    plot(ht, col=cols)
    abline(v=moda)
  }
  return(moda=moda)
}

ht <- hist(x)
moda.hist(ht)

#-----------------------------------------------------------------------------
# função para obter a moda da função densidade kernel

moda.dens <- function(dn, plotit=TRUE){
  ## dn: um objeto do uso da função density()
  ini <- dn$x[which.max(dn$y)]
  fx <- approxfun(dn$x, dn$y)
  op <- optim(c(ini), fx, method="BFGS", control=list(fnscale=-1))
  if(plotit==TRUE){
    plot(dn)
    abline(v=op$par, col=2)
  }
  return(moda=op$par)
}

dn <- density(x, kernel="triangular", bw=0.01)
moda.dens(dn)

dn <- density(x, kernel="gaussian", bw=0.05)
moda.dens(dn)

#-----------------------------------------------------------------------------
# gráfico

png("f027.png", 400, 300)
x <- rgamma(1000, 5, 5)
ht <- hist(x, freq=FALSE, col="gray86")
mo1 <- moda.hist(ht, plotit=FALSE)
abline(v=mo1, col=2, lwd=2)
dn <- density(x)
lines(dn, col=4, lwd=2)
mo2 <- moda.dens(dn, plotit=FALSE)
abline(v=mo2, col=3, lwd=2)
box()
dev.off()

#-----------------------------------------------------------------------------


Tags: , , , ,

Comments are closed.