Por que essa distribuição é uniforme?

12

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 pA=pB , ou seja, a probabilidade de conversão é igual nos dois grupos. Geramos dados artificiais usando um modelo binomial, por exemplo,

nABinomial(N,pA)

Em seguida, tentamos estimar pA,pB usando um modelo beta-binomial bayesiano para obter posteriores para cada taxa de conversão, por exemplo,

PABeta(1+nA,NnA+1)

Nossa estatística de teste é calculada calculando S=P(PA>PB|N,nA,nB) via monte carlo.

O que me surpreendeu foi que, se pA=pB , então SUniform(0,1) . Meu pensamento era que ele seria centrado em torno de 0,5 e até convergir para 0,5 conforme o tamanho da amostra, N , cresce.

Minha pergunta é: por que SUniform(0,1) quando pA=pB ?


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()
Cam.Davidson.Pilon
fonte
Observe que não pode ser exatamente uniforme, porque é uma variável discreta. Você está, portanto, perguntando sobre o comportamento assintótico. Além disso, para N pequeno (menos de 100 / min ( p , 1 - p ) , aproximadamente, comSN100/min(p,1p) ) a distribuição não é nem remotamente próxima do uniforme. p=pA=pB
whuber
@whuber S não é discreto, é uma probabilidade que pode cair entre 0 e 1. Além disso, mesmo para N baixo, estou observando um comportamento uniforme.
Cam.Davidson.Pilon
2
Devo estar entendendo mal sua configuração, então. Tanto quanto eu posso dizer, para qualquer valor dado de o valor de S é um número. Portanto, aceitando que N ,N,nA,nB,S e p B são fixadas para o momento (como eles estão em seu código), S é uma função de ( n A , n B ) . Mas a última, sendo realização de duas distribuições binomiais, pode atingir apenas um conjunto discreto de valores. Quando reproduzo seu código emN,pA,pBS(nA,nB)R, Obtenho histogramas decididamente não uniformes para pequeno . N
whuber
1
Embora seu tenha valores entre 0 e 1S01 , não confunda isso com não discreto: ele pode ter no máximo valores distintos (e na verdade tem menos que isso). Isso pode não estar perfeitamente claro para você, porque sua simulação gera estimativas de S em vez de seus valores corretos e as estimativas essencialmente têm uma distribuição contínua. N2S
whuber
1
@ whuber sim, você está correto, excelente observação. Eu ainda estou preso por que parece uniforme então.
Cam.Davidson.Pilon

Respostas:

11

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 ssimulationprimesimulationunderlyingsimulationprime .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).simulationprimeA = 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 usimulationunderlyingsimulationprime 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 gsimulationunderlyingBetaA>BetaBsimulationunderlying, 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:

a = b = 0.5
N = 10
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,1000)
P.show()
russellpierce
fonte
1
Portanto, há uma diferença entre o meu e o seu código. Eu amostro A e B em cada loop, você amostra uma vez e calcula S 5000 vezes.
Cam.Davidson.Pilon
1
A discrepância está nas suas chamadas para rbinom, que retornam um vetor. A chamada subseqüente para rbetadentro replicateé 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 ( ). ABNSIM = 10000?rbetaABreplicate
cardeal
1
aqui está a saída para os curiosos: imgur.com/ryvWbJO
Cam.Davidson.Pilon
1
As únicas coisas que sei que são potencialmente pertinentes em um nível conceitual são: a) a distribuição esperada dos resultados é simétrica, b) um tamanho de compartimento 1 é sempre uniforme, c) um tamanho de compartimento 2 para uma distribuição simétrica também sempre parecerá uniforme; d) o número de possíveis distribuições de amostragem que podem ser obtidas dos aumentos com N, e) os valores de S não podem se acumular apenas em 0 ou 1, porque beta é indefinido quando há 0 sucessos em ambos os grupos e f) as amostras são restritas entre 0 e 1.
russellpierce
1
Apenas como observação, podemos ver que as distâncias entre os centróides das distribuições de amostragem diminuem à medida que os centróides das distribuições de amostragem se afastam de 0,5 (provavelmente relacionado ao ponto f acima). Esse efeito tende a neutralizar a tendência das altas frequências de observações para os sucessos quase iguais mais comuns nos casos do grupo A e do grupo B. No entanto, fornecer uma solução matemática para saber por que isso é ou por que deve gerar distribuições normais para determinados tamanhos de lixeira não está nem perto do meu território.
russellpierce
16

NO(1/N)


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(1p)N. Consequently α=nA/N and β=nB/N have expectations of p and variance p(1p)/N.

  • As a Beta(nA+1,N+1nA) variate, PA has an expectation of (nA+1)/(N+2) and a variance of (nA+1)(N+1nA)/[(N+2)2(N+3)]. Approximating, we find that PA has an expectation of

    E(PA)=α+O(1/N)

    and a variance of

    Var(PA)=α(1α)/N+O(1/N2),

    with similar results for PB.

Let us therefore approximate the distributions of PA and PB with Normal(α,α(1α)/N) and Normal(β,β(1β)/N) distributions (where the second parameter designates the variance). The distribution of PAPB consequently is approximately Normal; to wit,

PAPBNormal(αβ,α(1α)+β(1β)N).

For very large N, the expression α(1α)+β(1β) will not vary appreciably from p(1p)+p(1p)=2p(1p) except with very low probability (another neglected O(1/N) term). Accordingly, letting Φ be the standard normal CDF,

Pr(PA>PB)=Pr(PAPB>0)Φ(αβ2p(1p)/N).

But since αβ has zero mean and variance 2p(1p)/N, Z=αβ2p(1p)/N is a standard Normal variate (at least approximately). Φ is its probability integral transform; Φ(Z) is uniform.

whuber
fonte
1
I'mn with you up until PAPBNormal... then you go off another direction that I didn't quite follow. Is Φ defined twice, once as the standard normal CDF and then as the probability integral transform? I'm hoping you can expand your description around these steps and relate them to the initial code/problem. Maybe loop back around and restate which specific parameters produce the uniform result.
russellpierce
1
@rpierce (1) The difference PAPB is approximately normal because PA and PB are independent and each is approximately normal. The mean is the difference of the means and the variance is the sum of the variances. (2) The probability integral transform is the CDF: it is the case for any random variable X with continuous distribution F, that F(X) is uniform.
whuber
1
Oh I got 1, it was the stuff after it where I got lost. This will be mindbogglingly dumb, but why is Pr(PA>PB) the same as the CDF?
russellpierce
1
@rpierce That follows rather directly from the definition, but there's a slight twist in which the symmetry of the Normal distribution is invoked. We're dealing with a Normal variate X=PAPB assumed to have an expectation of μ=αβ and variance σ2=2p(1p)/N. Standardizing X, it is natural to rewrite the probability as
Pr(X>0)=Pr((Xμ)/σ>(0μ)/σ)=1Φ(μ/σ)=Φ(μ/σ).
whuber
3
@whuber this is pretty amazing. You are a wonderful teacher. I appreciate both yours and rpierce's answer, I'll still give him credit as it did solve our problem, and you have shown why the behaviour occurs. Ty!
Cam.Davidson.Pilon