Estamos investigando las pruebas estadísticas bayesianas, y nos encontramos con un fenómeno extraño (al menos para mí).
Considere el siguiente caso: estamos interesados en medir qué población, A o B, tiene una tasa de conversión más alta. Para una comprobación de cordura, establecemos , es decir, la probabilidad de conversión es igual en ambos grupos. datos artificiales utilizando un modelo binomial, por ejemplo,n A ∼ Binomial ( N , p A )
Luego tratamos de estimar el usando un modelo beta-binomial bayesiano para obtener posteriores para cada tasa de conversión, por ejemplo,
Nuestra estadística de prueba se calcula calculando través de Monte Carlo.
Lo que me sorprendió fue que si , entonces . Mis pensamientos eran que sería centra alrededor de 0,5, e incluso convergen a 0,5 como el tamaño de la muestra, , crece.
Mi pregunta es, ¿por qué cuando ?
Aquí hay un código de Python para demostrar:
%pylab
from scipy.stats import beta
import numpy as np
import pylab as P
a = b = 0.5
N = 10000
samples = [] #collects the values of S
for i in range(5000):
assert a==b
A = np.random.binomial(N, a); B = np.random.binomial(N, b)
S = (beta.rvs(A+1, N-A+1, size=15000) > beta.rvs(B+1, N-B+1, size=15000)).mean()
samples.append(S)
P.hist(samples)
P.show()
fuente
R
Respuestas:
TL; DR: Las mezclas de distribuciones normales pueden verse uniformes cuando los tamaños de los contenedores son grandes.
Esta respuesta toma prestada del código de muestra de @ whuber (que primero pensé que era un error, pero en retrospectiva probablemente fue una pista).
Las proporciones subyacentes en la población son iguales:
a = b = 0.5
.Cada grupo, A y B tiene 10000 miembros:
N = 10000
.Vamos a realizar 5000 repeticiones de una simulación:
for i in range(5000):
.En realidad, lo que estamos haciendo es una de una . En cada una de las 5000 iteraciones haremos . s i m u l a t i o n u n d e r l y i n g s i m u l a t i o n p r i m esimulationprime simulationunderlying simulationprime simulationunderlying
En cada iteración del vamos a simular un número aleatorio de A y B que son 'éxitos' (también conocido como convertida) dadas las proporciones subyacentes iguales definidos anteriormente: . Nominalmente, esto producirá A = 5000 y B = 5000, pero A y B varían de una ejecución sim a otra y se distribuyen entre las 5000 ejecuciones de simulación de forma independiente y (aproximadamente) normalmente (volveremos a eso).simulationprime
A = np.random.binomial(N, a); B = np.random.binomial(N, b)
Pasemos ahora a para una única iteración de en la que A y B han tenido el mismo número de éxitos (como será el promedio). En cada iteración de , dados A y B, crearemos variantes aleatorias de la distribución beta para cada grupo. Luego los compararemos y descubriremos si , produciendo un VERDADERO o FALSO (1 o 0). Al final de una ejecución de , hemos completado 15000 iteraciones y tenemos 15000 valores VERDADERO / FALSO. El promedio de estos arrojará un solo valor de la distribución de muestreo (aproximadamente normal) de la proporción de s i m B e tsimulationunderlying s i m u l a t i o n u n d e r l y i n g B e t un a >simulationprime simulationunderlying simulatio n u n d e r l y i n g B e t a A > B e t a BBetaA>BetaB simulationunderlying BetaA>BetaB .
Excepto que ahora seleccionará los valores 5000 A y B. A y B rara vez serán exactamente iguales, pero las diferencias típicas en el número de éxitos de A y B se ven reducidas por el tamaño total de la muestra de A y B. Los As y Bs típicos producirán más tirones de su distribución de muestreo de proporciones de , pero los que están en los bordes de la distribución A / B también serán eliminados.B e t a A > B e t a Bsimulationprime BetaA>BetaB
Entonces, lo que en esencia detenemos muchas ejecuciones de simulación es una combinación de distribuciones de muestreo de para combinaciones de A y B (con más tirones de las distribuciones de muestreo hechas de los valores comunes de A y B que los valores poco comunes de A y B). Esto da como resultado mezclas de distribuciones normales-ish. Cuando los combina en un tamaño de contenedor pequeño (como es el valor predeterminado para la función de histograma que utilizó y se especificó directamente en su código original), termina con algo que parece una distribución uniforme.BetaA>BetaB
Considerar:
fuente
rbinom
, que devuelve un vector. La llamada posterior alrbeta
interiorreplicate
se vectoriza, por lo que el bucle interno (interno) está utilizando un y diferente para cada una de las 15000 variables aleatorias generadas (envolviendo las 5000 finales desde su ). Mira para más. Esto difiere del código de @ Cam con un solo y fijo utilizado en todas las 15000 llamadas de variable aleatoria para cada uno de los 5000 bucles de muestreo ( ). B A BNSIM = 10000
?rbeta
replicate
Para tener una idea de lo que está sucediendo, siéntase libre de hacer que muy grande y, al hacerlo, ignorar el comportamiento de y explotar los teoremas asintóticos que establecen que las distribuciones Beta y Binomial se vuelven aproximadamente normales. (Con algunos problemas, todo esto puede hacerse riguroso). Cuando hacemos esto, el resultado surge de una relación específica entre los diversos parámetros.O ( 1 / N )N O(1/N)
Debido a que planeamos usar aproximaciones normales, prestaremos atención a las expectativas y variaciones de las variables:
Como Binomial variables aleatorias, y tienen expectativas de y varianzas de . En consecuencia y tienen expectativas de y la varianza .n A n B p N p ( 1 - p ) N α = n A / N β = n B / N p p ( 1 - p ) / N(N,p) nA nB pN p(1−p)N α=nA/N β=nB/N p p(1−p)/N
Como una Beta , tiene una expectativa de y una varianza de . Aproximadamente, encontramos que tiene una expectativa de(nA+1,N+1−nA) PA (nA+1)/(N+2) (nA+1)(N+1−nA)/[(N+2)2(N+3)] PA
y una varianza de
con resultados similares para .PB
Por lo tanto, aproximaremos las distribuciones de y con Normal y Normal (donde el segundo parámetro designa la varianza ). La distribución de consecuencia es aproximadamente Normal; esto es,PA PB (α,α(1−α)/N) (β,β(1−β)/N) PA−PB
Para muy grande , la expresión no variará apreciablemente de excepto con muy baja probabilidad (otro término descuidado ). En consecuencia, dejando que sea el CDF normal estándar,N α(1−α)+β(1−β) p(1−p)+p(1−p)=2p(1−p) O(1/N) Φ
Pero dado que tiene media cero y varianza es una Normal estándar variante (al menos aproximadamente). es su probabilidad de transformación integral ; es uniforme .2 p ( 1 - p ) / N , Z = α - βα−β 2p(1−p)/N, ΦΦ(Z)Z=α−β2p(1−p)/N√ Φ Φ(Z)
fuente