¿Cuál es la relación entre el perfil de la probabilidad y los intervalos de confianza?

Question

Para hacer este cuadro me generaron muestras aleatorias de diferentes tamaños, desde una distribución normal con media=0 y dt=1. Los intervalos de confianza fueron calculados usando alfa atajos que van desde .001 a .999 (línea roja) con la t.función test (), el perfil de la probabilidad fue calculado utilizando el siguiente código que encontré en notas de la conferencia en línea (no puedo encontrar el enlace en el momento de Edición:Encontrado), esto se muestra por las líneas azules. Líneas verdes muestran la densidad normalizada mediante la R de la densidad() la función y los datos se muestran por el boxplots en la parte inferior de cada cuadro. En el derecho es una oruga trama del 95% intervalos de confianza (rojo) y 1/20 de max probabilidad de intervalos (azul).

R Código utilizado para el perfil de probabilidad:

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

Mi pregunta concreta es si existe una relación entre estos dos tipos de intervalos y por qué el intervalo de confianza parece ser más conservador para todos los casos, excepto cuando n=3. Comentarios/respuestas acerca de si mis cálculos son válidos (y una mejor manera de hacer esto) y la relación general entre estos dos tipos de intervalos.

R código:

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)

Answer 1

No voy a dar una respuesta completa (tengo un tiempo difícil tratando de entender lo que usted está haciendo exactamente), pero voy a tratar de aclarar cómo el perfil de probabilidad se construye. Puedo completar mi respuesta más tarde.

El pleno de la probabilidad de una muestra normal de tamaño $n$ es $$L(\mu, \sigma^2) = \left( \sigma^2 \right)^{-n/2} \exp\left( - \sum_i (x_i-\mu)^2/2\sigma^2 \right).$$

If $\mu$ is your parameter of interest, and $\sigma^2$ is a nuisance parameter, a solution to make inference only on $\mu$ es definir el perfil de probabilidad $$L_P(\mu) = L\left(\mu, \widehat{\sigma^2}(\mu) \right)$$ donde $\widehat{\sigma^2}(\mu)$ es el MLE para $\mu$ fijo: $$\widehat{\sigma^2}(\mu) = \text{argmax}_{\sigma^2} L(\mu, \sigma^2).$$

Se verifica que $$\widehat{\sigma^2}(\mu) = {1\over n} \sum_k (x_k - \mu)^2.$$

De ahí que el perfil de probabilidad es $$L_P(\mu) = \left( {1\over n} \sum_k (x_k - \mu)^2 \right)^{-n/2} \exp( -n/2 ).$$

Here is some R code to compute and plot the profile likelihood (I removed the constant term $\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")

Enlace con la probabilidad de que voy a tratar de resaltar el vínculo con la probabilidad de con el siguiente gráfico.

Primero definir la probabilidad:

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

A continuación, hacer un gráfico de contorno:

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

Y luego superponer la gráfica de $\widehat{\sigma^2}(\mu)$:

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

Los valores del perfil de probabilidad son los valores tomados por la probabilidad a lo largo de la red parábola.

Usted puede utilizar el perfil de la probabilidad como un univariante clásica de la probabilidad (cf @Prokofiev respuesta). Por ejemplo, el MLE $\hat\mu$ es el mismo.

Para su intervalo de confianza, los resultados difieren un poco debido a la curvatura de la función de $\widehat{\sigma^2}(\mu)$, pero en la medida en que se ocupan de un segmento corto de la misma, es casi lineal, y la diferencia va a ser muy pequeño.

También puede utilizar el perfil de probabilidad para construir la puntuación de las pruebas, por ejemplo.

Answer 2

En un marco general, el perfil de probabilidad de los intervalos aproximados de los intervalos de confianza. La prueba de este resultado es esencialmente el mismo como la demostración de que el cociente de probabilidad estadística (asintóticamente) distribuye aproximadamente como una $\chi^2_k$ distribución. La idea consiste en la inversión de la prueba de hipótesis obtenidos a partir de un cociente de probabilidad estadística.

Por ejemplo, $0.147$a nivel de perfil de la probabilidad de un intervalo, de un aproximado de confianza de $95\%$.

Estos son los resultados clásicos y, por tanto, me limitaré a ofrecer algunas referencias sobre esto:

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

El siguiente código R indica que, incluso para muestras pequeñas, los intervalos obtenidos con ambos métodos son similares (estoy re-uso de Elvis ejemplo):

Tenga en cuenta que usted tiene que utilizar la normalizado perfil de la probabilidad.

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

Si usamos un mayor tamaño de muestra, los intervalos de confianza son aún más cerca:

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

UN PUNTO IMPORTANTE:

Tenga en cuenta que para determinadas muestras de diferentes tipos de intervalos de confianza pueden diferir en términos de su longitud o ubicación, lo que realmente importa es su cobertura. En el largo plazo, todos ellos deben ofrecer la misma cobertura, independientemente de lo mucho que difieren de ejemplos específicos.