Calcular o intervalo de confiança para a média de uma distribuição beta

12

Considere uma distribuição beta para um determinado conjunto de classificações em [0,1]. Depois de calcular a média:

μ=αα+β

Existe uma maneira de fornecer um intervalo de confiança em torno dessa média?

dominic
fonte
1
dominic - você definiu a média da população . Um intervalo de confiança seria baseado em alguma estimativa dessa média. Que estatística de amostra você está usando?
Glen_b -Reinstate Monica
Glen_b - Olá, estou usando um conjunto de classificações normalizadas (de um produto) no intervalo [0,1]. O que estou procurando é uma estimativa de um intervalo em torno da média (para um determinado nível de confiança), por exemplo: média + - 0,02
dominic
2
±
Glen_b - obrigado pela sua paciência. Vou usar o MLE
dominic
2
nμμ

Respostas:

22

Embora existam métodos específicos para calcular intervalos de confiança para os parâmetros em uma distribuição beta, descreverei alguns métodos gerais, que podem ser usados ​​para (quase) todos os tipos de distribuições , incluindo a distribuição beta, e são facilmente implementados em R .

Intervalos de confiança da probabilidade do perfil

Vamos começar com a estimativa de probabilidade máxima com intervalos de confiança de probabilidade de perfil correspondentes. Primeiro, precisamos de alguns dados de amostra:

# Sample size
n = 10

# Parameters of the beta distribution
alpha = 10
beta = 1.4

# Simulate some data
set.seed(1)
x = rbeta(n, alpha, beta)

# Note that the distribution is not symmetrical
curve(dbeta(x,alpha,beta))

Função densidade de probabilidade para a distribuição beta.

A média real / teórica é

> alpha/(alpha+beta)
0.877193

Agora precisamos criar uma função para calcular a função de probabilidade de log negativa para uma amostra da distribuição beta, com a média como um dos parâmetros. Podemos usar a dbeta()função, mas como isso não usa uma parametrização envolvendo a média, temos que expressar seus parâmetros ( α e β ) como uma função da média e algum outro parâmetro (como o desvio padrão):

# Negative log likelihood for the beta distribution
nloglikbeta = function(mu, sig) {
  alpha = mu^2*(1-mu)/sig^2-mu
  beta = alpha*(1/mu-1)
  -sum(dbeta(x, alpha, beta, log=TRUE))
}

Para encontrar a estimativa de probabilidade máxima, podemos usar a mle()função na stats4biblioteca:

library(stats4)
est = mle(nloglikbeta, start=list(mu=mean(x), sig=sd(x)))

Apenas ignore os avisos por enquanto. Eles são causados ​​pelos algoritmos de otimização que tentam valores inválidos para os parâmetros, fornecendo valores negativos para α e / ou β . (Para evitar o aviso, você pode adicionar um lowerargumento e alterar a otimização methodusada.)

Agora, temos estimativas e intervalos de confiança para nossos dois parâmetros:

> est
Call:
mle(minuslogl = nloglikbeta, start = list(mu = mean(x), sig = sd(x)))

Coefficients:
        mu        sig 
0.87304148 0.07129112

> confint(est)
Profiling...
         2.5 %    97.5 %
mu  0.81336555 0.9120350
sig 0.04679421 0.1276783

Observe que, como esperado, os intervalos de confiança não são simétricos:

par(mfrow=c(1,2))
plot(profile(est)) # Profile likelihood plot

Gráfico de probabilidade de perfil para a distribuição beta.

(As segundas linhas magentas externas mostram o intervalo de confiança de 95%.)

Observe também que, mesmo com apenas 10 observações, obtemos estimativas muito boas (um intervalo de confiança estreito).

Como alternativa mle(), você pode usar a fitdistr()função do MASSpacote. Isso também calcula o estimador de probabilidade máxima e tem a vantagem de que você só precisa fornecer a densidade, não a probabilidade negativa do log, mas não fornece intervalos de confiança da probabilidade do perfil, apenas intervalos de confiança assintóticos (simétricos).

Uma opção melhor é mle2()(e funções relacionadas) do bbmlepacote, que é um pouco mais flexível e poderoso do que mle()e oferece gráficos um pouco mais agradáveis.

Intervalos de confiança de inicialização

Outra opção é usar o bootstrap. É extremamente fácil de usar no R, e você nem precisa fornecer uma função de densidade:

> library(simpleboot)
> x.boot = one.boot(x, mean, R=10^4)
> hist(x.boot)                # Looks good
> boot.ci(x.boot, type="bca") # Confidence interval
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = x.boot, type = "bca")

Intervals : 
Level       BCa          
95%   ( 0.8246,  0.9132 )  
Calculations and Intervals on Original Scale

O bootstrap tem a vantagem adicional de funcionar, mesmo que seus dados não venham de uma distribuição beta.

Intervalos de confiança assintóticos

Para intervalos de confiança médios, não vamos esquecer os bons e velhos intervalos de confiança assintóticos baseados no teorema do limite central (e na distribuição t ). Desde que tenhamos um tamanho de amostra grande (para que o CLT se aplique e a distribuição da média da amostra seja aproximadamente normal) ou grandes valores de α e β (para que a distribuição beta em si seja aproximadamente normal), ela funcionará bem. Aqui não temos nenhum, mas o intervalo de confiança ainda não é tão ruim:

> t.test(x)$conf.int
[1] 0.8190565 0.9268349

Para valores ligeiramente maiores de n (e não valores extremos demais para os dois parâmetros), o intervalo de confiança assintótica funciona extremamente bem.

Karl Ove Hufthammer
fonte
Obrigado Karl. Pergunta rápida: como você determinou seu alfa e beta? Usei a variância e a média da amostra para obter alfa e beta, mas acho que posso ter confundido a média da amostra com a média da população, por isso não tenho certeza de ter feito isso da maneira certa ... veja o comentário de Glen_b acima .
dominic
Para determinar α e β como funções da média e desvio padrão, acabei de inverter as funções dos desvios médios e padrão como funções de α e β (mas tenho certeza de que você também pode pesquisar na rede).
Karl Ove Hufthammer
α,β
0

Confira a regressão beta. Uma boa introdução a como fazê-lo usando R pode ser encontrada aqui:

http://cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf

Outra maneira (realmente fácil) de construir um intervalo de confiança seria usar uma abordagem não paramétrica de boostrap. A Wikipedia tem boas informações:

http://en.wikipedia.org/wiki/Bootstrapping_%28statistics%29

Também um bom vídeo aqui:

http://www.youtube.com/watch?v=ZCXg64l9R_4

Rasmus Bååth
fonte