Incapaz de reproducir el documento de resultados (tamaño de muestra)

Question

Estoy interesado en utilizar los resultados de este trabajo (Modificado exacto del tamaño de la muestra para una proporción binomial con especial énfasis en el diagnóstico resto de estimación de parámetros por Geoffrey Fosgate) para calcular los tamaños de muestra. Existe una implementación del algoritmo utilizado para calcular el tamaño de la muestra en R binomSamSize paquete. He comprobado el código que se usa para calcular esto y parece seguir el papel exactamente, pero los resultados que se muestran en el papel de la sección de resultados no parecen ser reproducible. El tamaño de la muestra llegar a ser demasiado pequeño.

El papel requiere una suscripción para ser visto en su totalidad, pero aquí es un fragmento de código para mostrar cómo el tamaño de la muestra se calcula:

Aquí están los resultados que desea reproducir:

Este es el código que estoy utilizando (levantado de la binomSamSize paquete):

ciss.midp <- function(p0, d, alpha, nMax = 1e+06){
  pi.L <- p0 - d
  pi.U <- p0 + d
  if (pi.L < 0) 
    stop("p0 - d is below zero!")
  if (pi.U > 1) 
    stop("p0 + d is above one!")
  n <- floor(max(1/p0, 1/(1 - p0)))
  done <- FALSE
  while (!done & (n < nMax)) {
    n <- n + 1
    x <- round(p0 * n)
    lhs2 <- 1/2 * dbinom(x, size = n, prob = pi.L) + 1/2 * 
      dbinom(x, size = n, prob = pi.U) + 
      pbinom(x, size = n, prob = pi.L, lower.tail = FALSE) + 
      pbinom(x - 1, size = n, prob = pi.U)
    if (!is.na(lhs2)) {
      done <- (lhs2 < alpha)
    }
  }
  return(n)
}

Lo que me acaban de llegar a producir un círculo columna es este:

> sapply(seq(0.5, 0.9, 0.05), function(i) ciss.midp(p0=i, d=0.1, alpha=0.1))
[1] 68 67 65 61 57 50 42 34 23

PS: Si hay algo más que pueda para hacerlo más fácil de responder, por favor hágamelo saber en los comentarios.

Answer 1

Tal vez el autor del aplicación - a pesar de la descripción - se inicia con un gran n y las obras de la espalda y se detiene en la primera vez que estamos por encima de alfa, por lo tanto devolución de los últimos n donde la suma estaba por debajo de la alfa? Al menos un la implantación de un algo que sería consistente con los resultados de la de papel (código de prueba a continuación).

Parecen existir más de una n, donde un interruptor de la suma que se por encima de alfa de n y estar por debajo del alfa de n+1... sin Embargo, yo creo que es mejor ir desde pequeñas n y de seguridad.

##Go the other way...not my 1st choice though...
ciss.midp2 <- function (p0, d, alpha, nMax = 1e+05) {
    pi.L <- p0 - d
    pi.U <- p0 + d
    if (pi.L < 0)
        stop("p0 - d is below zero!")
    if (pi.U > 1)
        stop("p0 + d is above one!")
    n <- nMax
    done <- FALSE
    while (!done & (n > 0)) {
        n <- n - 1
        x <- round(p0 * n)
        lhs2 <- 1/2 * dbinom(x, size = n, prob = pi.L) +
          pbinom(x, size = n, prob = pi.L, lower.tail = FALSE) +
          (1 - pbinom(x-1, size = n, prob = pi.U, lower.tail=FALSE)) +
          1/2*dbinom(x, size = n, prob = pi.U)

        if (!is.na(lhs2)) {
            done <- (lhs2 > alpha)
        }
    }
    return(n+1)
}

library("binomSamSize")
p.grid <- seq(0.5, 0.9, 0.05)
sapply(p.grid, function(i) ciss.midp(p0=i, d=0.1, alpha=0.1))
[1] 68 67 65 61 57 50 42 34 23
sapply(p.grid, function(i) ciss.midp2(p0=i, d=0.1, alpha=0.1)
[1] 68 67 65 61 57 52 45 39 29