Estimar el tamaño de la muestra si se miden dos propiedades y se conocen sus probabilidades

Question

Lo siento si esto es un ejercicio estándar, pero no sé cómo buscarlo.

Tengo una cantidad infinita de muestras y puedo medir dos propiedades $A$ y $B$ para cada muestra. Posibilidad de medir $A$ es $P(A)$ y para $B$ es $P(B)$ ambas propiedades son independientes. Tomo N muestras y en $a$ de los que mido la propiedad $A$ pero la propiedad $B$ se encuentra $b$ tiempos.

Si $P(A)$ , $P(B)$ , $a$ y $b$ son conocidos, ¿cómo puedo conseguir $N$ y ¿cuál es el error/desviación de esta cantidad?

También puedo formular la pregunta como un problema con algunos números: Tengo una cantidad infinita de cajas, en cada caja hay un caramelo envuelto en papel. El 34% de los caramelos son de chocolate y el 10% de los papeles son azules.

He abierto algunas cajas y he visto 16 papeles azules y me he comido 45 chocolates. ¿Cuántas cajas he abierto?

¿Qué cambia si $P(A)$ y $P(B)$ no son precisos, sino que se miden experimentalmente a partir de $N_A$ y $N_B$ ¿medidas en su lugar? Eso significa que - antes de mi medidas mi hermano abrió algunos $N_A=522399$ cajas y ahí es donde el $P(A)=34%$ vino de. Es daltónico, pero mi hermana lo vio abriendo algunos $N_B=92587$ cajas y noté que el 10% de los papeles eran azules. ¿Algún cambio si hay más de dos propiedades?

Answer 1

Abordaremos este problema utilizando estimación de máxima verosimilitud para estimar $N$ .

Primero, escribamos la probabilidad de ver $x$ chocolates y $y$ papeles azules, entregados $N$ cajas fueron abiertas, la probabilidad de ver un chocolate es $p_x$ y la probabilidad de ver un papel azul es $p_y$ . Esto no es más que el producto de dos distribuciones binomiales con el mismo $N$ :

$P(x,y) = \frac{N!}{x!(N-x)!} p_x^x(1-p_x)^{N-x}\frac{N!}{y!(N-y)!}p_y^y(1-p_y)^{N-y}$

Resolvemos el $N$ que nos da la mayor probabilidad de ver las observaciones que realmente vimos, en este caso 45 chocolates y 16 papeles azules. La mejor manera de hacerlo es maximizar el logaritmo de $P(x,y)$ por razones de estabilidad y para evitar números realmente grandes o pequeños.

Aquí hay algo de código R para hacer esto de una manera bastante bruta:

log.pxy <- function(N, px, py, x, y)
{
  dbinom(x, N, px, log=TRUE) + dbinom(y, N, py, log=TRUE)
}

results <- rep(-Inf,200)
for (N in 50:200) results[N] <- log.pxy(N, 0.34, 0.1, 45, 16)
which.max(results)
[1] 137

Nuestra estimación de máxima probabilidad de $N$ es 137.

En cuanto a la obtención de un intervalo de confianza, podemos utilizar el hecho de que -2 * el logaritmo observado de la función de verosimilitud ( $P(x,y)$ visto como una función de $N$ en lugar de $(x,y)$ ) se distribuye asintóticamente $\chi^2(1)$ . Encontramos el $N$ para los que la probabilidad logarítmica está "demasiado lejos" por debajo del máximo, utilizando el $\chi^2(1)$ como guía, y construyendo un intervalo de confianza del 95%:

> min(which(2*results > max(2*results)-qchisq(0.95,1)))
[1] 111
> max(which(2*results > max(2*results)-qchisq(0.95,1)))
[1] 168

Por lo tanto, un intervalo de confianza del 95% sería [110, 169] - tenemos que ampliar el intervalo en 1 en cada extremo para llegar "fuera" del qchisq Rango del 95%.

En cuanto a sus otras preguntas: Si hay más de dos propiedades, puedes ampliar la metodología de la solución de la manera obvia y seguirá funcionando.

La situación más compleja es cuando $p_x$ y $p_y$ son a su vez estimaciones basadas en muestras. Yo estimaría conjuntamente $p_x$ , $p_y$ y $N$ en ese caso, lo que puede hacerse con un procedimiento anidado que itere sobre $N$ en un bucle externo, y, en una optimización interna, estima $p_x$ y $p_y$ dados todos los datos y $N$ . (Como se verá, esto es bastante sencillo; dado $N$ Es sólo la estimación habitual de $p_x$ y $p_y$ .) A continuación, encontramos el $N$ que maximiza la probabilidad logarítmica como antes.

La función de probabilidad es más compleja. Vamos a denotar la otra información con $N_a$ y $N_b$ tamaños de muestra y valores observados $x_a$ y $y_b$ . Tenemos:

$L(N,p_x,p_y) = {{N}\choose{x}}p_x^x(1-p_x)^{N-x} {{N_a}\choose{x_a}}p_x^{x_a}(1-p_x)^{N_a-x_a} \dots$

donde el $\dots$ nos ahorra la redacción de la $p_y$ parte. Obviamente, podemos combinar algunos términos, pero esta forma hace que sea un poco más fácil ver lo que está pasando.

Ahora el código R. Asumiré, para concretar, que $N_a=200$ y $x_a=68$ , dando la estimación puntual de $p_x=0.34$ y $N_b=100$ , $y_b=10$ , dando la estimación puntual de $p_y=0.1$ .

log.ll <- function(px, py, N, x, y, Na, Nb, xa, yb) {
  dbinom(x, N, px, log=TRUE) + dbinom(y, N, py, log=TRUE) +
      dbinom(xa, Na, px, log=TRUE) + dbinom(yb, Nb, py, log=TRUE)
}

x = 45
y = 16
Na = 200
Nb = 100
xa = 68
yb = 10

log.ll.N <- rep(-Inf,200)
for (N in 51:200) {
  px.hat <- (x+xa)/(N+Na)
  py.hat <- (y+yb)/(N+Nb)

  log.ll.N[N] <- log.ll(px.hat, py.hat, N, x, y, Na, Nb, xa, yb)
}

Y, por las respuestas:

> which.max(log.ll.N)
[1] 135
> min(which(2*log.ll.N > max(2*log.ll.N)-qchisq(0.95,1)))
[1] 103
> max(which(2*log.ll.N > max(2*log.ll.N)-qchisq(0.95,1)))
[1] 180
>

Para una estimación puntual ligeramente diferente de 135 para $N$ y un intervalo de confianza más amplio de 102 - 181, como corresponde a nuestra nueva falta de precisión sobre $p_x$ y $p_y$ . Podemos recuperar nuestras nuevas estimaciones de $p_x$ y $p_y$ basado en nuestra muestra combinada:

> N <- 135
> (x+xa)/(N+Na)
[1] 0.3373134
> (y+yb)/(N+Nb)
[1] 0.1106383

También debo señalar que nuestro intervalo de confianza se basa en el perfil log probabilidad no la probabilidad logarítmica, pero sigue siendo un intervalo de confianza perfectamente válido.

Answer 2

Esto me recuerda al marcar y recapturar estrategia para estimar el tamaño de una población. Esto está modelado por la distribución hipergeométrica, pero en su caso cada propiedad es una binomial.

Para una propiedad se escribiría la probabilidad de observar $a$ como

$${{N}\choose{a}} P(A)^a(1-P(A))^{N-a}.$$

Como las propiedades son independientes se puede escribir toda la probabilidad como el producto

$$L(N) = {{N}\choose{a}} P(A)^a(1-P(A))^{N-a} {{N}\choose{b}} P(B)^k(1-P(B))^{N-b},$$

que hay que maximizar para $N$ . Yo resolvería esto usando R porque tomando el logaritmo y diferenciando con respecto a $N$ introduce la función digamma (la derivada de la función log-Gamma), que no simplifica las cosas.

Aquí hay un R script que daría la estimación.

p <- .34  # P(A)
q <- .1   # P(B)
a <- 45
b <- 16
# Minimum possible value for N.
N <- max(a,b)
old_loglik <- -Inf
while (TRUE) {
   loglik <- dbinom(a, size=N, prob=p, log=TRUE) +
       dbinom(b, size=N, prob=q, log=TRUE)
   if (old_loglik > loglik) {
      # Passed the optimum.
      N <- N-1
      break
   }
   else {
      old_loglik <- loglik
   }
   N <- N+1
}
print (N)
# Answer is 137.

Si $P(A)$ y $P(B)$ no se conocen y hay que estimar su valor, entonces el problema inicial carece de sentido porque se necesita saber $N$ para estimarlos (que alguien me corrija si conoce un estimador alternativo).

Si se estimaron en una primera muestra de tamaño $n$ y se extrae una nueva muestra de tamaño $N$ entonces el problema es más complicado porque hay que escribir la probabilidad de todo y maximizarla, esta vez con 3 parámetros $(p, q, N)$ .