Estamos investigando testes estatísticos bayesianos e deparamos com um fenômeno estranho (para mim pelo menos).
Considere o seguinte caso: estamos interessados em medir qual população, A ou B, tem uma taxa de conversão mais alta. Para uma verificação de sanidade, definimos , ou seja, a probabilidade de conversão é igual nos dois grupos. Geramos dados artificiais usando um modelo binomial, por exemplo,
Em seguida, tentamos estimar usando um modelo beta-binomial bayesiano para obter posteriores para cada taxa de conversão, por exemplo,
Nossa estatística de teste é calculada calculando via monte carlo.
O que me surpreendeu foi que, se , então . Meu pensamento era que ele seria centrado em torno de 0,5 e até convergir para 0,5 conforme o tamanho da amostra, , cresce.
Minha pergunta é: por que quando ?
Aqui está um código Python para demonstrar:
%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()
fonte
R
, Obtenho histogramas decididamente não uniformes para pequeno .Respostas:
TL; DR: misturas de distribuições normais podem parecer uniformes quando os tamanhos de lixeira são grandes.
Essa resposta é emprestada do código de exemplo do @ whuber (que eu pensei que era primeiro um erro, mas em retrospecto provavelmente era uma dica).
As proporções subjacentes na população são iguais:
a = b = 0.5
.Cada grupo, A e B tem 10000 membros:
N = 10000
.Nós estamos indo para conduzir 5000 repetições de uma simulação:
for i in range(5000):
.Na realidade, o que está a fazer é um de um s i m u l a t i o n u n d e r L y i n g . Em cada um dos 5000 iterações s i m u l a t i o n p r i m e faremos ssimulationprime simulationunderlying simulationprime .simulationunderlying
Em cada iteração de que vai simular um número aleatório de A e B que são 'sucessos' (AKA convertido) tendo em conta as proporções subjacentes iguais definidas anteriormente:. Nominalmente, isso produzirá A = 5000 e B = 5000, mas A e B variam de execução sim para execução sim e são distribuídos pelas corridas da simulação 5000 de forma independente e (aproximadamente) normalmente (voltaremos a isso).simulationprime
A = np.random.binomial(N, a); B = np.random.binomial(N, b)
Vamos agora passo através para uma única iteração de s i m u l a t i o n p r i m e em que A e B têm obtido em igual número de sucessos (como será a média do caso). Em cada iteração de s i m u l a t i o n usimulationunderlying simulationprime iremos, dado A e B, criar variates aleatórias da distribuição beta para cada grupo. Então, vamos compará-los e descobrir se B e t um A > B e t um B , produzindo um VERDADEIRO ou FALSO (1 ou 0). No final de uma corrida desimulatio n u n d e r L y i n gsimulationunderlying BetaA>BetaB simulationunderlying , concluímos 15000 iterações e temos 15000 valores VERDADEIRO / FALSO. A média destes irá produzir um valor único de amostragem a distribuição (aproximadamente normal) da proporção de .BetaA>BetaB
Só que agora vai seleccionar 5000 valores de A e B. A e B raramente serão exatamente iguais, mas as diferenças típicas no número de sucessos de A e B são diminuídas pelo tamanho total da amostra de A e B. As e Bs típicas renderão mais puxões de sua distribuição amostral de proporções de B e t asimulationprime , mas aqueles nas bordas da distribuição A / B também serão puxados.BetaA>BetaB
Então, o que, em essência, nós encostar várias execuções sim é uma combinação de amostragem distribuições de para combinações de A e B (com mais puxa a partir das distribuições de amostragem feitas a partir dos valores comuns de um e B que os valores incomuns de A e B). Isso resulta em misturas de distribuições normais. Quando você os combina em um tamanho de compartimento pequeno (como é o padrão para a função de histograma usada e especificada diretamente no código original), você acaba com algo que se parece com uma distribuição uniforme.BetaA>BetaB
Considerar:
fonte
rbinom
, que retornam um vetor. A chamada subseqüente pararbeta
dentroreplicate
é vetorizada, de modo que o loop interno (interno) está usando A e B diferentes para cada uma das 15000 variáveis aleatórias geradas (agrupando os 5000 finais desde o seu ). Veja para mais. Isso difere do código do @ Cam, com um único A e B fixo usado em todas as 15000 chamadas de variável aleatória para cada um dos 5000 loops de amostragem ( ).NSIM = 10000
?rbeta
replicate
Como planejamos usar aproximações normais, prestaremos atenção às expectativas e variações das variáveis:
As Binomial(N,p) variates, nA and nB have expectations of pN and variances of p(1−p)N . Consequently α=nA/N and β=nB/N have expectations of p and variance p(1−p)/N .
As a Beta(nA+1,N+1−nA) variate, PA has an expectation of (nA+1)/(N+2) and a variance of (nA+1)(N+1−nA)/[(N+2)2(N+3)] . Approximating, we find that PA has an expectation of
and a variance of
with similar results forPB .
Let us therefore approximate the distributions ofPA and PB with Normal(α,α(1−α)/N) and Normal(β,β(1−β)/N) distributions (where the second parameter designates the variance). The distribution of PA−PB consequently is approximately Normal; to wit,
For very largeN , the expression α(1−α)+β(1−β) will not vary appreciably from p(1−p)+p(1−p)=2p(1−p) except with very low probability (another neglected O(1/N) term). Accordingly, letting Φ be the standard normal CDF,
But sinceα−β has zero mean and variance 2p(1−p)/N, Z=α−β2p(1−p)/N√ is a standard Normal variate (at least approximately). Φ is its probability integral transform; Φ(Z) is uniform.
fonte