Jaki jest związek między prawdopodobieństwem profilu a przedziałami ufności?

18

Aby zrobić ten wykres, wygenerowałem losowe próbki o różnej wielkości z rozkładu normalnego ze średnią = 0 i sd = 1. Przedziały ufności zostały następnie obliczone przy użyciu wartości odcięcia alfa w zakresie od 0,001 do .999 (czerwona linia) za pomocą funkcji t.test (), prawdopodobieństwo profilu zostało obliczone przy użyciu kodu poniżej, który znalazłem w notatkach z wykładów umieszczonych w linii (mogę t znajdź link w tej chwili Edytuj: Znaleziono ), jest to pokazane przez niebieskie linie. Zielone linie pokazują znormalizowaną gęstość za pomocą funkcji gęstości R (), a dane są pokazane przez wykresy prostokątne u dołu każdego wykresu. Po prawej stronie znajduje się wykres gąsienicowy z 95% przedziałami ufności (czerwony) i 1/20 maksymalnych przedziałów prawdopodobieństwa (niebieski).

Kod R użyty dla prawdopodobieństwa profilu:

  #mn=mean(dat)
  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )

wprowadź opis zdjęcia tutaj

Moje konkretne pytanie dotyczy tego, czy istnieje znany związek między tymi dwoma typami przedziałów i dlaczego przedział ufności wydaje się bardziej zachowawczy dla wszystkich przypadków, z wyjątkiem sytuacji, gdy n = 3. Pożądane są również komentarze / odpowiedzi na temat tego, czy moje obliczenia są prawidłowe (i lepszy sposób na zrobienie tego) oraz ogólny związek między tymi dwoma typami przedziałów.

Kod R:

samp.size=c(3,4,5,10,20,1000)
cnt2<-1
ints=matrix(nrow=length(samp.size),ncol=4)
layout(matrix(c(1,2,7,3,4,7,5,6,7),nrow=3,ncol=3, byrow=T))
par(mar=c(5.1,4.1,4.1,4.1))
for(j in samp.size){


  #set.seed(200)
  dat<-rnorm(j,0,1)
  vals<-seq(.001,.999, by=.001)
  cis<-matrix(nrow=length(vals),ncol=3)
  cnt<-1
  for(ci in vals){
    x<-t.test(dat,conf.level=ci)$conf.int[1:2]
    cis[cnt,]<-cbind(ci,x[1],x[2])
    cnt<-cnt+1
  }


  mn=mean(dat)
  n=length(dat)
  high<-max(c(dat,cis[970,3]), na.rm=T)
  low<-min(c(dat,cis[970,2]), na.rm=T)
  #high<-max(abs(c(dat,cis[970,2],cis[970,3])), na.rm=T)
  #low<--high


  muVals <- seq(low,high, length = 1000)
  likVals <- sapply(muVals,
                    function(mu){
                      (sum((dat - mu)^2) /
                         sum((dat - mn)^2)) ^ (-n/2)
                    }
  )


  plot(muVals, likVals, type = "l", lwd=3, col="Blue", xlim=c(low,high),
       ylim=c(-.1,1), ylab="Likelihood/Alpha", xlab="Values",
       main=c(paste("n=",n), 
              "True Mean=0 True sd=1", 
              paste("Sample Mean=", round(mn,2), "Sample sd=", round(sd(dat),2)))
  )
  axis(side=4,at=seq(0,1,length=6),
       labels=round(seq(0,max(density(dat)$y),length=6),2))
  mtext(4, text="Density", line=2.2,cex=.8)

  lines(density(dat)$x,density(dat)$y/max(density(dat)$y), lwd=2, col="Green")
  lines(range(muVals[likVals>1/20]), c(1/20,1/20), col="Blue", lwd=4)
  lines(cis[,2],1-cis[,1], lwd=3, col="Red")
  lines(cis[,3],1-cis[,1], lwd=3, col="Red")
  lines(cis[which(round(cis[,1],3)==.95),2:3],rep(.05,2), 
        lty=3, lwd=4, col="Red")
  abline(v=mn, lty=2, lwd=2)
  #abline(h=.05, lty=3, lwd=4, col="Red")
  abline(h=0, lty=1, lwd=3)
  abline(v=0, lty=3, lwd=1)

  boxplot(dat,at=-.1,add=T, horizontal=T, boxwex=.1, col="Green")
  stripchart(dat,at=-.1,add=T, pch=16, cex=1.1)

  legend("topleft", legend=c("Likelihood"," Confidence Interval", "Sample Density"),
         col=c("Blue","Red", "Green"), lwd=3,bty="n")

  ints[cnt2,]<-cbind(range(muVals[likVals>1/20])[1],range(muVals[likVals>1/20])[2],
                     cis[which(round(cis[,1],3)==.95),2],cis[which(round(cis[,1],3)==.95),3])
  cnt2<-cnt2+1
}
par(mar=c(5.1,4.1,4.1,2.1))


plot(0,0, type="n", ylim=c(1,nrow(ints)+.5), xlim=c(min(ints),max(ints)), 
     yaxt="n", ylab="Sample Size", xlab="Values")
for(i in 1:nrow(ints)){
  segments(ints[i,1],i+.2,ints[i,2],i+.2, lwd=3, col="Blue")
  segments(ints[i,3],i+.3,ints[i,4],i+.3, lwd=3, col="Red")
}
axis(side=2, at=seq(1.25,nrow(ints)+.25,by=1), samp.size)
Kolba
źródło
W notatkach z wykładów mnjest literówka mu, a nie mean(dat). Jak powiedziałem w komentarzach do drugiego pytania , powinno to być jasne na stronie 23 z definicjami.
Elvis
@Elvis Nie sądzę. mn jest zdefiniowany na stronie 18 notatek.
Flask
Próbowałem wyjaśnić pojęcie prawdopodobieństwa profilu. Czy możesz skomentować nieco dalej, co robisz w powyższym kodzie?
Elvis
3
@Elvis Ani też nie rozumiem. Przedział ufności oparty na prawdopodobieństwie profilu powinien być konstruowany za pomocą percentyli, które nigdzie się nie pojawiają. χ2
Stéphane Laurent,
1
@ StéphaneLaurent nie jestem pewien oryginalny kod jest zapewnienie przedziały ufności. Raczej 1/20 maksymalne przedziały prawdopodobieństwa. Wierzę, że nazwa przedziałów ufności na moim wykresie to przedziały ufności typu „wald”, a czerwone linie na wykresach to „krzywe zaufania” opisane na tej stronie wikipedii
Flask

Odpowiedzi:

10

Nie udzielę pełnej odpowiedzi (trudno mi zrozumieć, co dokładnie robisz), ale postaram się wyjaśnić, w jaki sposób budowane jest prawdopodobieństwo profilu. Mogę uzupełnić odpowiedź później.

Pełne się prawdopodobieństwo prawidłowej próbki o wymiarach jest L ( μ , σ 2 ) = ( σ 2 ) - n / 2 exp ( - Σ I ( x i - μ ) 2 / 2 σ 2 ) .n

L(μ,σ2)=(σ2)n/2exp(i(xiμ)2/2σ2).

μσ2μ

LP(μ)=L(μ,σ2^(μ))
σ2^(μ)μ
σ2^(μ)=argmaxσ2L(μ,σ2).

σ2^(μ)=1nk(xkμ)2.

LP(μ)=(1nk(xkμ)2)n/2exp(n/2).

exp(n/2)

> data(sleep)
> difference <- sleep$extra[11:20]-sleep$extra[1:10]
> Lp <- function(mu, x) {n <- length(x); mean( (x-mu)**2 )**(-n/2) }
> mu <- seq(0,3, length=501)
> plot(mu, sapply(mu, Lp, x = difference), type="l")

prawdopodobieństwo profilu

Link z prawdopodobieństwem Spróbuję podświetlić link z prawdopodobieństwem na poniższym wykresie.

Najpierw określ prawdopodobieństwo:

L <- function(mu,s2,x) {n <- length(x); s2**(-n/2)*exp( -sum((x-mu)**2)/2/s2 )}

Następnie wykonaj wykres konturowy:

sigma <- seq(0.5,4, length=501)
mu <- seq(0,3, length=501)

z <- matrix( nrow=length(mu), ncol=length(sigma))
for(i in 1:length(mu))
  for(j in 1:length(sigma))
    z[i,j] <- L(mu[i], sigma[j], difference)

# shorter version
# z <- outer(mu, sigma, Vectorize(function(a,b) L(a,b,difference)))

contour(mu, sigma, z, levels=c(1e-10,1e-6,2e-5,1e-4,2e-4,4e-4,6e-4,8e-4,1e-3,1.2e-3,1.4e-3))

σ2^(μ)

hats2mu <- sapply(mu, function(mu0) mean( (difference-mu0)**2 ))
lines(mu, hats2mu, col="red", lwd=2)

wykres konturowy L.

Wartości prawdopodobieństwa profilu to wartości przyjęte przez prawdopodobieństwo wzdłuż czerwonej paraboli.

μ^

σ2^(μ)

Możesz również użyć prawdopodobieństwa profilu do zbudowania testów wyników, na przykład.

Elvis
źródło
mu w kodzie jest sekwencją wartości od niskiej do wysokiej, prawdopodobieństwo dla każdej z tych wartości jest dzielone przez prawdopodobieństwo dla średniej próbki (mn). Jest to znormalizowane prawdopodobieństwo.
Flask
Myślę, że to to samo, ale nie znormalizowane. Czy możesz umieścić go w kodzie R lub w inny sposób wykreślić funkcję dla niektórych danych, abyśmy mogli porównać?
Flask
Oto jest Na początku myślałem, że mnto literówka, teraz myślę, że kod R jest nieprawidłowy. Sprawdzę to jutro dwukrotnie - żyję już późno.
Elvis
Możesz mieć rację. Nie rozumiem, jak kod udaje się go znormalizować. Och, rozumiem, „normalizacja” dzieli się tylko przez maksimum?
Elvis
1
Myślę, że powinno to ułatwić dostrzeżenie, kiedy iloraz prawdopodobieństwa jest mniejszy niż pewien próg (np. 1/20 maks.) Przy jakiejś hipotezie zerowej (np. Zero).
Flask
7

χk2

0.14795%

Są to klasyczne wyniki, dlatego przedstawię kilka odniesień na ten temat:

http://www.jstor.org/stable/2347496

http://www.stata-journal.com/sjpdf.html?articlenum=st0132

http://www.unc.edu/courses/2010fall/ecol/563/001/docs/lectures/lecture11.htm

http://en.wikipedia.org/wiki/Likelihood-ratio_test

http://en.wikipedia.org/wiki/Likelihood_function#Profile_likelihood

Poniższy kod R pokazuje, że nawet w przypadku małych próbek przedziały uzyskane przy obu podejściach są podobne (ponownie używam przykładu Elvisa):

Pamiętaj, że musisz użyć znormalizowanego prawdopodobieństwa profilu.

data(sleep)
x <- sleep$extra[11:20]-sleep$extra[1:10]
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(0,3, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(0.5,1.5))$root,uniroot(Rpt,c(1.51,3))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

Jeśli użyjemy większej próby, przedziały ufności będą jeszcze bliższe:

set.seed(123)
x <- rnorm(100)
n <- length(x)
Rp <- function(mu) {mean( (x-mean(x))^2 )^(n/2)/mean( (x-mu)^2 )^(n/2) }
Rp(mean(x))

mu <- seq(-0.5,0.5, length=501)
plot(mu, sapply(mu, Rp), type="l")


Rpt<- function(mu) Rp(mu)-0.147 # Just an instrumental function

# Likelihood-confidence interval of 95% level

c(uniroot(Rpt,c(-0.4,0))$root,uniroot(Rpt,c(0,0.4))$root)

# t confidence interval

t.test(x,conf.level=0.95)$conf.int

WAŻNY PUNKT:

Należy pamiętać, że w przypadku konkretnych próbek różne rodzaje przedziałów ufności mogą się różnić pod względem długości lub lokalizacji, tak naprawdę liczy się ich zasięg. Na dłuższą metę wszystkie powinny zapewniać taki sam zasięg, niezależnie od tego, jak bardzo różnią się dla poszczególnych próbek.

Prokofiew
źródło
@Prokoflev, jeśli istnieje jakaś prosta zależność między przedziałami ufności obliczonymi za pomocą funkcji R t.test () a tymi obliczonymi na podstawie powyższego kodu funkcji prawdopodobieństwa, możesz to opublikować. Szczególnie interesuje mnie przypadek n = 3. Niestety mam niewielkie doświadczenie w matematyce, więc wiele artykułów prowadzi mnie do króliczej nory, szukając nazw symboli i ich reprezentacji itp., Gdy kilka wierszy kodu (najłatwiejszy jest R) może mi to wyjaśnić.
Flask
@Flask Czy jesteś zainteresowany uzyskaniem przedziałów ufności dla parametrów rozkładu normalnego lub bardziej ogólnych ram?
Prokofiew
@Prokoflev specjalnie dla średniej rozkładu normalnego, jak pokazano w moim przykładzie w pytaniu. Szczególnie zastanawiam się, dlaczego przedziały ufności są bardziej zachowawcze, z wyjątkiem przypadku n = 3.
Flask
95%
1
Zaczynam wierzyć, że powinienem pomnożyć przedziały prawdopodobieństwa przez jakieś kwantyle rozkładu normalnego lub chisquare, aby uzyskać odpowiedni przedział ufności.
Flask
1

χ2)normzaljazmire

  1. Prawdopodobieństwo dziennika profilu jest w przybliżeniu kwadratowe
  2. Istnieje transformacja parametru, która sprawia, że ​​prawdopodobieństwo logarytmu profilu jest w przybliżeniu kwadratowe.

Kwadrat jest ważny, ponieważ określa normalny rozkład w skali logarytmicznej. Im bardziej kwadratowy, tym lepsze przybliżenie i wynikające CI. Twój wybór 1/20 wartości granicznej dla przedziałów prawdopodobieństwa jest równoważny ponad 95% CI w granicy asymptotycznej, dlatego właśnie niebieskie przedziały są na ogół dłuższe niż czerwone.

Jest jeszcze jeden problem z prawdopodobieństwem profilu, który wymaga uwagi. Jeśli masz wiele zmiennych, których profilujesz, to jeśli liczba punktów danych na wymiar jest niska, prawdopodobieństwo profilu może być bardzo stronnicze i optymistyczne. W celu zmniejszenia tego obciążenia wykorzystuje się zatem marginalne, warunkowe i zmodyfikowane prawdopodobieństwa profilu.

Tak więc odpowiedź na twoje pytanie brzmi TAK ... związek jest asymptotyczną normalnością większości estymatorów maksymalnego prawdopodobieństwa, co przejawia się w rozkładzie chi-kwadrat współczynnika prawdopodobieństwa.


źródło
Jeśli masz wiele zmiennych, których profilujesz, to jeśli liczba punktów danych na wymiar jest niska, prawdopodobieństwo profilu może być bardzo stronnicze i optymistyczne ” Optymistyczny w porównaniu z czym?
Flask
@Flask Przez optymizm rozumiem, że będzie on zbyt wąski, aby zapewnić nominalne prawdopodobieństwo pokrycia przy traktowaniu go jako przedziału ufności.
Rozumiem, dziękuję, ale w moim konkretnym przypadku jest to rzeczywiście pesymistyczne? Nie jestem pewien, czy mówimy o przedziałach prawdopodobieństwa, czy przedziałach ufności wynikających z prawdopodobieństw.
Flask
@Flask Myślę, że interwały wydają się pesymistyczne, ponieważ porównujesz 1/20 przedział wiarygodności (5% względnego prawdopodobieństwa) z 95% CI. Jak twierdzą inni tutaj, naprawdę chciałbyś porównać to z 15% względnym przedziałem prawdopodobieństwa, aby jabłka od jabłek ... przynajmniej asymptotycznie. Twój przedział prawdopodobieństwa na obecnym etapie rozważa więcej opcji jako wiarygodne.
Mam szczegółowe rzeczywistego problemu pragnę zastosować co uczę się tutaj . Martwię się, że w przypadku, gdy rozkład próbkowania jest nieznany (ale prawdopodobnie nie normalny) i skomplikowany, twoje dwa wymagania mogą nie zostać spełnione. Jednak obliczone przeze mnie prawdopodobieństwo profilu wydaje się normalne i rozsądne. Czy to dlatego, że rozkład próbkowania średniej powinien być normalnie rozłożony?
Flask