0

Necesito obtener los percentiles 5 % y 95% en la salida de la instrucción bugs en R2openBUGS, para poder hacer un contraste de hipótesis unilateral, con un intervalo de credibilidad bayesiano del 95%. Les adjunto el código en r y openBUGS. El problema es que no encuentro la instrucción dentro de Bugs en R2openBUGS para poder indicar o setear más percentiles o cuantiles. ¿Pueden ustedes por favor ayudarme en esta tarea? Muchas gracias por la ayuda.

        library("R2OpenBUGS") 
        library("coda") 
        library("lattice")

        n4 = 10000  
        n5 = 1000
        n6 = 1
        n7 = 1
        N= 10
        r = 8
        a = 1
        b = 1

        mod<-"mbp2022.txt"

        data <- list( N = N,
              r = r,
              a = a,
              b = b          
        )
        sink(mod)
        cat("
       
        model{

       pi1 ~ dbeta(a,b)    
       r ~ dbin(pi1, N)

      }

      ",fill=TRUE) 
      sink()

      modelfile <- mod


      vars <-list("pi1")

      inits <-function(){list(pi1=0.5)} 

 

      out <-bugs(
  
      model.file=file.path(modelfile),
  
      data=data,
  
      inits=inits,
  
      parameters.to.save=vars, 
   
      n.iter=n4, 
      n.burnin=n5, 
      n.chains=n6,
      n.thin=n7,
      summary.only = F,
      debug=F, 
      DIC=F, 
      working.directory = getwd(),
  
      )
      cat("\n")  
      print (out, digits = 03)

1 respuesta 1

0

Para encontrar los percentiles se puede utilizar la función quantile pero depende de qué data quieras. No conozco en detalle la salida de openBugs pero supón que la salida out$sims.matrix es una matriz con iteraciones en filas y variables simuladas en columnas y la segunda columna es la variable buscada:

quantile(out$sims.matrix[,2],c(0.05,0.95))

devuelve los cuantiles 5% y 95%

7
  • La solución aportada es buena, funciona solo cuando tengo más de 1 variable, para cada fila la salida es correcta. Ahora, por ejemplo la salida de la variable out (print (out, digits = 03) mean sd 2.5% 25% 50% 75% 97.5% pi1 0.748 0.121 0.481 0.674 0.765 0.841 0.928 di 0.248 0.121 -0.019 0.174 0.265 0.341 0.428 pp 0.033 0.180 0.000 0.000 0.000 0.000 1.000 Lo que necesito es el intervalo %5 y 95% de las tres filas o variables al mismo tiempo. He intentado manejando índices pero no lo logro.
    – toto28
    Commented el 7 nov. 2022 a las 13:22
  • Es que lo que me estás mostrando ya son cuantiles calculados, y no incluyen el 5% y el 95%. ¿algún parámetro para incluir esos cuantiles en BUGS? ¿Esto que imprimes es un resumen o el detalle de la variable? ¿que te muestra la salida de dput(out)?
    – Ric
    Commented el 7 nov. 2022 a las 13:30
  • Exacto. No encuentro ningún parámetro en BUGS para añadir esos %, y la salida de dput(out) no la puedo poner aquí porque es demasiado larga
    – toto28
    Commented el 7 nov. 2022 a las 14:08
  • Si out es un vector (p.e. unclass(out)) basta con quantile(unclass(out),c(0.05,0.95)), pero no lo se.
    – Ric
    Commented el 7 nov. 2022 a las 14:17
  • quantile(unclass(out),c(0.05,0.95)) me da error Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) : 'x' must be atomic y ` out$sims.list $pi1 [1] 0.5793 0.6720 0.8380 0.5790 es la salida de out$sims.list`
    – toto28
    Commented el 7 nov. 2022 a las 14:40

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.