Rozważ rozkład beta dla danego zestawu ocen w [0,1]. Po obliczeniu średniej:

Czy istnieje sposób na zapewnienie przedziału ufności wokół tego środka?

Chociaż istnieją określone metody obliczania przedziałów ufności dla parametrów w rozkładzie beta, opiszę kilka ogólnych metod, które można zastosować do (prawie) wszystkich rodzajów rozkładów , w tym rozkładu beta, i które można łatwo wdrożyć w R .

Przedziały ufności prawdopodobieństwa profilu

Zacznijmy od oszacowania maksymalnego prawdopodobieństwa z odpowiednimi przedziałami ufności prawdopodobieństwa profilu. Najpierw potrzebujemy przykładowych danych:

# 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))

Rzeczywista / teoretyczna średnia to

> alpha/(alpha+beta)
0.877193

Teraz musimy stworzyć funkcję do obliczania funkcji prawdopodobieństwa ujemnego dziennika dla próbki z rozkładu beta, ze średnią jako jednym z parametrów. Możemy użyć tej dbeta()funkcji, ale ponieważ nie używa ona parametryzacji obejmującej średnią, musimy wyrazić jej parametry ( α i β ) jako funkcję średniej i kilku innych parametrów (takich jak odchylenie standardowe):

# 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))
}

Aby znaleźć oszacowanie maksymalnego prawdopodobieństwa, możemy użyć mle()funkcji w stats4bibliotece:

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

Na razie zignoruj ​​ostrzeżenia. Są one spowodowane przez algorytmy optymalizujące próbujące nieprawidłowe wartości parametrów, dające wartości ujemne dla α i / lub β . (Aby uniknąć ostrzeżenia, możesz dodać lowerargument i zmienić zastosowaną optymalizację method).

Teraz mamy zarówno szacunki, jak i przedziały ufności dla naszych dwóch parametrów:

> 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

Należy pamiętać, że zgodnie z oczekiwaniami przedziały ufności nie są symetryczne:

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

(Druga zewnętrzna magenta pokazuje 95% przedział ufności.)

Zauważ też, że nawet przy zaledwie 10 obserwacjach otrzymujemy bardzo dobre szacunki (wąski przedział ufności).

Alternatywnie mle()możesz użyć fitdistr()funkcji z MASSpakietu. To również oblicza estymator maksymalnego prawdopodobieństwa i ma tę zaletę, że wystarczy podać gęstość, a nie ujemne prawdopodobieństwo dziennika, ale nie daje przedziałów ufności profilu, tylko asymptotyczne (symetryczne) przedziały ufności.

Lepszą opcją jest mle2()(i powiązane funkcje) z bbmlepakietu, który jest nieco bardziej elastyczny i wydajny niż mle()i daje nieco ładniejsze wykresy.

Przedziały ufności Bootstrap

Inną opcją jest użycie bootstrapu. Jest bardzo łatwy w użyciu w R i nie musisz nawet podawać funkcji gęstości:

> 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

Bootstrap ma tę dodatkową zaletę, że działa, nawet jeśli dane nie pochodzą z wersji beta.

Asymptotyczne przedziały ufności

W przypadku przedziałów ufności dla średniej, nie zapominajmy o starych dobrych asymptotycznych przedziałach ufności opartych na centralnym twierdzeniu granicznym (i rozkładzie t ). Tak długo, jak mamy duży rozmiar próbki (więc obowiązuje CLT i rozkład średniej próbki jest w przybliżeniu normalny) lub duże wartości zarówno α, jak i β (tak, że sam rozkład beta jest w przybliżeniu normalny), działa dobrze. Tutaj nie mamy żadnego, ale przedział ufności wciąż nie jest taki zły:

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

W przypadku nieznacznie dużych wartości n (i niezbyt ekstremalnych wartości dwóch parametrów) asymptotyczny przedział ufności działa wyjątkowo dobrze.

Sprawdź regresję beta. Dobre wprowadzenie do tego, jak to zrobić za pomocą R, można znaleźć tutaj:

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

Innym (naprawdę łatwym) sposobem konstruowania przedziału ufności byłoby zastosowanie nieparametrycznego podejścia przypominającego. Wikipedia ma dobre informacje:

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

Również fajne wideo tutaj:

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