Desvio padrão de observações em bin

24

Eu tenho um conjunto de dados de observações de amostra, armazenadas como contagens dentro de intervalos. por exemplo:

min/max  count
40/44    1
45/49    2
50/54    3
55/59    4
70/74    1

Agora, encontrar uma estimativa da média disso é bastante simples. Basta usar a média (ou mediana) de cada compartimento de intervalo como observação e a contagem como peso e encontrar a média ponderada:

x¯=1i=1Nwii=1Nwixi

Para o meu caso de teste, isso me dá 53,82.

Minha pergunta agora é: qual é o método correto para encontrar o desvio padrão (ou variação)?

Através da minha pesquisa, encontrei várias respostas, mas não tenho certeza qual, se houver, é realmente apropriado para o meu conjunto de dados. Consegui encontrar a fórmula a seguir em outra pergunta aqui e em um documento NIST aleatório .

s2=i=1Nwi(xix¯)2(M1)Mi=1Nwi

O que fornece um desvio padrão de 8,35 para o meu caso de teste. No entanto, o artigo da Wikipedia sobre médias ponderadas fornece a fórmula:

s2=i=1Nwi(i=1Nwi)2i=1Nwi2i=1Nwi(xix¯)2

e

s2=1 1(Eu=1 1NWEu)-1 1Eu=1 1NWEu(xEu-x¯)2

Que fornecem desvios padrão de 8,66 e 7,83, respectivamente, para o meu caso de teste.

Atualizar

Agradecemos a @whuber que sugeriu analisar as correções de Sheppard e seus comentários úteis relacionados a elas. Infelizmente, estou tendo dificuldades para entender os recursos que posso encontrar sobre isso (e não consigo encontrar bons exemplos). Para recapitular, entendo que o seguinte é uma estimativa tendenciosa da variação:

s2=1 1Eu=1 1NWEuEu=1 1NWEu(xEu-x¯)2

Também entendo que a maioria das correções padrão para o viés são para amostras aleatórias diretas de uma distribuição normal. Portanto, vejo dois problemas em potencial para mim:

  1. Essas são amostras aleatórias em bin (que, tenho certeza, é onde as correções de Sheppard entram).
  2. Não se sabe se os dados são ou não para uma distribuição normal (portanto, suponho que não, o que, tenho certeza, invalida as correções de Sheppard).

Então, minha pergunta atualizada é; Qual é o método apropriado para lidar com o viés imposto pela fórmula "simples" de desvio / variação padrão ponderada em uma distribuição não normal? Mais especificamente no que diz respeito aos dados em bin.

Nota: estou usando os seguintes termos:

  • s2 é a variação ponderada
  • N é o número de observações. (ou seja, o número de posições)
  • M é o número de pesos diferentes de zero. (ou seja, o número de posições com contagens)
  • WEu são os pesos (ou seja, as contagens)
  • xEu são as observações. (ou seja, o bin significa)
  • x¯ é a média ponderada.
chezy525
fonte
5
Google "Correções de Sheppard" para as soluções padrão para esse problema.
whuber
@whuber, receio que meu google-foo esteja falhando comigo ... não estou descobrindo muito sobre como usar as correções de Sheppard. Tanto quanto posso dizer, é uma correção para a natureza dos dados em binários, e no meu caso de teste seria usado como , em que é o tamanho das caixas (no meu caso de teste, 4). Isso está correto? De qualquer forma, o que estou descobrindo ainda não parece me ajudar com a computação . cs2*s2-c212cs2
Chezy525
O segundo hit na minha pesquisa no Google fornece uma fórmula explícita (equação 9).
whuber
@whuber, já faz alguns meses, e eu tentei ler o documento que você vinculou algumas vezes. Acho que ainda estou perdendo alguma coisa, mas o melhor que descobri é que a equação final que listei está correta como estimador imparcial. Isto está certo?
chezy525
As correções de Sheppard não assumem normalidade.
Glen_b -Reinstala Monica

Respostas:

18

Esta resposta apresenta duas soluções: as correções de Sheppard e uma estimativa de máxima verossimilhança. Ambos concordam com uma estimativa do desvio padrão: para o primeiro e para o segundo (quando ajustado para ser comparável ao estimador "imparcial" usual).7,697,707.69


Correções de Sheppard

"Correções de Sheppard" são fórmulas que ajustam momentos calculados a partir de dados em bin (como estes) em que

  • presume-se que os dados sejam governados por uma distribuição suportada em um intervalo finito[uma,b]

  • esse intervalo é dividido sequencialmente em compartimentos iguais de largura comum que é relativamente pequeno (nenhum compartimento contém uma grande proporção de todos os dados)h

  • a distribuição tem uma função de densidade contínua.

Eles são derivados da fórmula da soma de Euler-Maclaurin, que aproxima integrais em termos de combinações lineares de valores do integrando em pontos regularmente espaçados e, portanto, geralmente aplicáveis ​​(e não apenas às distribuições normais).

Embora estritamente falando, uma distribuição Normal não seja suportada em um intervalo finito, é muito aproximada. Essencialmente, toda a sua probabilidade está contida em sete desvios-padrão da média. Portanto, as correções de Sheppard são aplicáveis ​​aos dados assumidos como provenientes de uma distribuição Normal.

As duas primeiras correções de Sheppard são

  1. Use a média dos dados em bin para a média dos dados (ou seja, nenhuma correção é necessária para a média).

  2. Subtraia 2/12 da variação dos dados em bin para obter a variação (aproximada) dos dados.h2/12

De onde vem 2/12? Isso é igual à variação de uma variável uniforme distribuída ao longo de um intervalo de comprimento . Intuitivamente, a correção de Sheppard para o segundo momento sugere que o armazenamento em cache dos dados - efetivamente substituindo-os pelo ponto médio de cada compartimento - parece acrescentar um valor aproximadamente uniformemente distribuído, variando entre e , de onde inflaciona a variação por .h - h / 2 h / 2h2/12h-h/2h/2h2/12

Vamos fazer os cálculos. Eu uso Rpara ilustrá-los, começando especificando as contagens e as caixas:

counts <- c(1,2,3,4,1)
bin.lower <- c(40, 45, 50, 55, 70)
bin.upper <- c(45, 50, 55, 60, 75)

A fórmula adequada a ser usada para as contagens vem da replicação das larguras dos compartimentos pelas quantidades fornecidas pelas contagens; isto é, os dados binados são equivalentes a

42.5, 47.5, 47.5, 52.5, 52.5, 57.5, 57.5, 57.5, 57.5, 72.5

Seu número, média e variância podem ser calculados diretamente sem a necessidade de expandir os dados desta maneira: quando um compartimento tem o ponto médio e uma contagem de , então sua contribuição para a soma dos quadrados é . Isso leva à segunda das fórmulas da Wikipedia citadas na pergunta.k k x 2xkkx2

bin.mid <- (bin.upper + bin.lower)/2
n <- sum(counts)
mu <- sum(bin.mid * counts) / n
sigma2 <- (sum(bin.mid^2 * counts) - n * mu^2) / (n-1)

A média ( mu) é (sem necessidade de correção) e a variação ( ) é . (Sua raiz quadrada é conforme indicado na pergunta.) Como a largura do compartimento comum é , subtraímos da variação e pegamos sua raiz quadrada, obtendo para o desvio padrão.1195/2254,32sigma2675/1161,367,83h=5h2/12=25/122.08675/11-52/127,70


Estimativas de máxima verossimilhança

Um método alternativo é aplicar uma estimativa de probabilidade máxima. Quando a distribuição subjacente assumida tem uma função de distribuição (dependendo dos parâmetros a serem estimados) e o compartimento contém valores de um conjunto de valores independentes e idênticos distribuídos de , então o contribuição (aditiva) para a probabilidade de log desse bin é θ ( x 0 , x 1 ] k F θFθθ(x0 0,x1 1]kFθ

registroEu=1 1k(Fθ(x1 1)-Fθ(x0 0))=kregistro(Fθ(x1 1)-Fθ(x0 0))

(consulte MLE / Probabilidade de intervalo lognormalmente distribuído ).

A soma de todos os compartimentos fornece a probabilidade de log para o conjunto de dados. Como sempre, encontramos uma estimativa que minimiza . Isso requer otimização numérica e é acelerada fornecendo bons valores iniciais para . O código a seguir faz o trabalho para uma distribuição Normal:Λ(θ)θ^-Λ(θ)θR

sigma <- sqrt(sigma2) # Crude starting estimate for the SD
likelihood.log <- function(theta, counts, bin.lower, bin.upper) {
  mu <- theta[1]; sigma <- theta[2]
  -sum(sapply(1:length(counts), function(i) {
    counts[i] * 
      log(pnorm(bin.upper[i], mu, sigma) - pnorm(bin.lower[i], mu, sigma))
  }))
}
coefficients <- optim(c(mu, sigma), function(theta) 
  likelihood.log(theta, counts, bin.lower, bin.upper))$par

Os coeficientes resultantes são .(μ^,σ^)=(54,32,7.33)

Lembre-se, no entanto, que para distribuições normais, a estimativa de probabilidade máxima de (quando os dados são fornecidos exatamente e não empilhados) é o SD da população dos dados, e não a estimativa "corrigida de viés" mais convencional na qual a variação é multiplicada por . Vamos então (para comparação) corrigir o MLE de , encontrando . Isso se compara favoravelmente com o resultado da correção de Sheppard, que foi de .σn/(n-1 1)σn/(n-1 1)σ^=11/10×7.33=7.697,70

Verificando as premissas

Para visualizar esses resultados, podemos plotar a densidade normal ajustada sobre um histograma:

hist(unlist(mapply(function(x,y) rep(x,y), bin.mid, counts)),
     breaks = breaks, xlab="Values", main="Data and Normal Fit")
curve(dnorm(x, coefficients[1], coefficients[2]), 
      from=min(bin.lower), to=max(bin.upper), 
      add=TRUE, col="Blue", lwd=2)

Figura

Para alguns, isso pode não parecer um bom ajuste. No entanto, como o conjunto de dados é pequeno (apenas valores), podem ocorrer desvios surpreendentemente grandes entre a distribuição das observações e a verdadeira distribuição subjacente.11

Vamos verificar formalmente a suposição (feita pelo MLE) de que os dados são governados por uma distribuição Normal. Um teste aproximado de qualidade do ajuste pode ser obtido a partir de um : os parâmetros estimados indicam a quantidade esperada de dados em cada compartimento; a estatística compara as contagens observadas com as contadas esperadas. Aqui está um teste em :χ2χ2R

breaks <- sort(unique(c(bin.lower, bin.upper)))
fit <- mapply(function(l, u) exp(-likelihood.log(coefficients, 1, l, u)),
              c(-Inf, breaks), c(breaks, Inf))
observed <- sapply(breaks[-length(breaks)], function(x) sum((counts)[bin.lower <= x])) -
  sapply(breaks[-1], function(x) sum((counts)[bin.upper < x]))
chisq.test(c(0, observed, 0), p=fit, simulate.p.value=TRUE)

A saída é

Chi-squared test for given probabilities with simulated p-value (based on 2000 replicates)

data:  c(0, observed, 0) 
X-squared = 7.9581, df = NA, p-value = 0.2449

O software realizou um teste de permutação (necessário porque a estatística do teste não segue exatamente uma distribuição qui-quadrado: veja minha análise em Como entender os graus de liberdade ). Seu valor-p de , que não é pequeno, mostra muito pouca evidência de desvio da normalidade: temos motivos para confiar nos resultados de máxima probabilidade.0,245

whuber
fonte