Sejam variáveis aleatórias independentes assumindo valores ou com probabilidade 0,5 cada. Considere a soma . Desejo limitar a probabilidade . O melhor limite que tenho agora é onde c é uma constante universal. Isso é obtido através do limite inferior da probabilidade Pr (| x_1 + \ dots + x_n | <\ sqrt {t}) e Pr (| y_1 + \ dots + y_n | <\ sqrt {t}) pela aplicação de limites simples de Chernoff. Posso esperar obter algo significativamente melhor do que esse limite? Para iniciantes, posso pelo menos obter . Se eu conseguir caudas sub-gaussianas, isso provavelmente seria o melhor, mas podemos esperar isso (acho que não, mas não consigo pensar em um argumento)?
probability
random-variable
bernoulli-distribution
user1189053
fonte
fonte
Respostas:
A relação algébrica
exibe como o produto de duas somas independentes. Como e são variáveis independentes de Bernoulli , é uma variável Binomial que foi duplicado e alterado. Portanto, sua média é e sua variação é . Da mesma forma tem uma média de e variância de . Vamos padronizá-los agora, definindoS (xi+1)/2 (yj+1)/2 (1/2) X=∑ai=1xi (a,1/2) 0 a Y=∑bj=1yj 0 b
de onde
Com um alto (e quantificável) grau de precisão, à medida que cresce, aproxima-se da distribuição normal padrão. Portanto, vamos aproximar como vezes o produto de duas normais normais.a Xa S ab−−√
O próximo passo é perceber que
é um múltiplo da diferença dos quadrados de variáveis independentes normal padrão e . A distribuição de pode ser calculada analiticamente ( invertendo a função característica ): seu pdf é proporcional à função Bessel da ordem zero, . Porque esta função tem caudas exponenciais, podemos concluir imediatamente que a grande e e fixa , não há melhor aproximação para que dado na questão.U V Zab K0(|z|)/π a b t Pra,b(S>t)
Resta algum espaço para melhorias quando um (pelo menos) de e não é grande ou em pontos na cauda do perto de . Cálculos diretos da distribuição de mostram uma redução gradual das probabilidades da cauda em pontos muito maiores que , aproximadamente além de . Esses gráficos log-lineares do CDF de para vários valores de (dados nos títulos) (variando aproximadamente os mesmos valores que , distinguidos pela cor em cada gráfico) mostram o que está acontecendo. Para referência, o gráfico do limitea b S ±ab S ab−−√ abmax(a,b)−−−−−−−−−−√ S a b a K0 a distribuição é mostrada em preto. (Como é simétrico em torno de , , basta observar a cauda negativa.)S 0 Pr(S>t)=Pr(−S<−t)
À medida que cresce, o CDF se aproxima da linha de referência.b
Caracterizar e quantificar essa curvatura exigiria uma análise mais fina da aproximação Normal às variáveis binomiais.
A qualidade da aproximação da função de Bessel se torna mais clara nessas partes ampliadas (no canto superior direito de cada gráfico). Nós já estamos muito distantes. Embora a escala vertical logarítmica pode esconder diferenças substanciais, claramente pelo tempo atingiu a aproximação é bom para .a 500 |S|<ab√
Código R para calcular a distribuição deS
O seguinte levará alguns segundos para ser executado. (Ele calcula vários milhões de probabilidades para 36 combinações de e .) Em máquinas mais lentas, omitir as maiores um ou dois valores de e e aumentar o limite de trama inferior de para cerca de .a b 10−300 10−160
a
b
fonte
1/2 (1 + y BesselK[0,-y] StruveL[-1, y] - y BesselK[1,-y] StruveL[0, y])
. Seria interessante ver como: (a) o limite do OP executa, e (b) sua aproximação Normal, para o caso que observamos acima, ou seja, derivado usando a solução discreta exata do pmf.Comentário: editei o título na tentativa de refletir melhor que tipo de RVs são considerados na pergunta. Qualquer pessoa pode reeditar.
Motivação: Eu acho que não há necessidade de aceitar um limite superior, se pudermos derivar a distribuição de. ( ATUALIZAÇÃO : Não podemos ver os comentários e a resposta de Whuber).|Sab|
Denote . É fácil verificar que 's têm a mesma distribuição que o ' s e o 's. A função geradora de momento éZk=XiYj,k=1,...,ab Z X Y
Além disso, os são, a princípio, independentes em pares: A variável (os índices podem ser qualquer um), tem suporte com probabilidades correspondentes . Sua função de geração de momento éZ W=Z1+Z2 {−2,0,2} {1/4,1/2,1/4}
suspeitar que a independência total é válida, como segue (é óbvio para os mais sábios?): Para esta parte, indique . Então, pela regra da cadeiaZij=XiYj
Pela independência entre pares, temos . Considere . e são condicionais independentes em portanto, temos a segunda igualdade pela independência entre pares. Mas isso implica queP[Z12∣Z11]=P[Z12]
P[Z13,Z12∣Z11] Z13 Z12 Z11
Etc (eu acho). ( ATUALIZAÇÃO : Acho errado . Independência provavelmente vale para qualquer trigêmeo, mas não para todo o grupo. Então, o que se segue é apenas a derivação da distribuição de uma simples caminhada aleatória, e não uma resposta correta para a pergunta - veja Wolfies e Respostas de Whuber).
Se a independência total realmente se mantiver, temos a tarefa de derivar a distribuição de uma soma de
que parece uma simples caminhada aleatória , embora sem a clara interpretação deste último como uma sequência.
Se o suporte de será o número inteiro par em incluindo zero, enquanto se o suporte de será o número inteiro ímpar em , sem zero.ab=even S [−ab,...,ab] ab=odd S [−ab,...,ab]
Tratamos o caso de . Indique como o número de assumem o valor . Então o suporte de pode ser escrito . Para qualquer , obtemos um valor único para . Além disso, devido a probabilidades simétricas e independência (ou apenas trocabilidade?), Todas as realizações conjuntas possíveis das variáveis são equivalentes. Então, contamos e descobrimos que a função de massa de probabilidade de é,ab=odd
m Z −1 S S∈{ab−2m;m∈Z+∪{0};m≤ab} m S Z {Z1=z1,...,Zab=zab} S
Definindo , e número ímpar por construção, e o elemento típico do suporte de , temoss≡ab−2m S
Movendo-se para, como se , a distribuição de é simétrica em torno de zero, sem alocar a massa de probabilidade para zero e, portanto, a distribuição deé obtido "dobrando" o gráfico de densidade em torno do eixo vertical, dobrando essencialmente as probabilidades de valores positivos,|S| ab=odd S |S|
Então a função de distribuição é
Portanto, para qualquer real , , obtemos a probabilidade requeridat 1≤t<ab
Observe que a indicação garante que a soma será executada apenas até os valores incluídos no suporte de- por exemplo, se estabelecermos , ainda será executado até , uma vez que é obrigado a ser estranho, além de ser um inteiro.i=odd |S| t=10.5 i 9
fonte
Não é uma resposta, mas um comentário sobre a interessante resposta de Alecos, que é muito longa para caber em uma caixa de comentários.
Sejam variáveis aleatórias independentes do Rademacher e sejam variáveis aleatórias independentes do Rademacher. Alecos observa que:(X1,...,Xa) (Y1,...,Yb)
"... parece um simples passeio aleatório ". Se fosse como um simples passeio aleatório, a distribuição de seria simétrica 'unimodal em forma de sino' em torno de 0.S
Para ilustrar que é não um simples passeio aleatório, aqui está uma rápida comparação Monte Carlo de:
Claramente, não é uma simples caminhada aleatória; Observe também que S não é distribuído em todos os números pares (ou ímpares).S
Monte Carlo
Aqui é o código (em Mathematica ) utilizado para gerar uma única iteração da soma , dado e :S a b
Em seguida, 500.000 tais caminhos, dizer quando e , pode ser gerada com:a=5 b=7
O domínio de apoio para esta combinação de e é:a b
fonte
a
eb
menos de 1000, pelo menos) comorademacher[a_] := Transpose[{Range[-a, a, 2], Array[Binomial[a, #] &, a + 1, 0] /2^a}]; s[a_, b_] := {#[[1, 1]], Total[#[[;; , 2]]]} & /@ GatherBy[Flatten[Outer[Times, rademacher[a], rademacher[b], 1], 1], First]; ListLogPlot[s[5, 7]]
Experimente, digamoss[100,211]
,.WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]]
... é duas vezes mais rápido que aOuter
abordagem. Curioso para saber qual código você está usando? [Ambas as abordagens podem, é claro, ser feita mais rapidamente utilizandoParallelTable
, etc]sum[n_, a_, b_] := Block[{w, p}, w[x_] := Array[Binomial[x, #] &, x + 1, 0] /2^x; p[x_] := RandomChoice[w[x] -> Range[-x, x, 2], n]; p[a] p[b]]
. Então horaTally[sum[500000, 5, 7]]
. ParaR
aficianodos, a seguir faz a mesma coisa e leva apenas 50% mais do que o Mathematica :s <- function(n, a, b) (2 * rbinom(n, a, 1/2) - a)*(2 * rbinom(n, b, 1/2) - b); system.time(x <- table(s(5*10^5, 5, 7))); plot(log(x), col="#00000020")
.