Soma dos produtos das variáveis ​​aleatórias Rademacher

9

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 obterx1xa,y1yb+11S=i,jxi×yjP(|S|>t)2ectmax(a,b)cPr(|x1++xn|<t)Pr(|y1++yn|<t)ectab . 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)?

user1189053
fonte
Você já pensou em aplicar o Chernoff vinculado diretamente a S ? Você pode fazer algo com
E[exp(λS]=E[λijXiYj]=E[λ(iXi)(jYj)]
Dilip Sarwate
Há uma melhoria óbvia no seu limite para t>ab , pois a probabilidade deve ser zero. Parece-me que é uma cauda "sub-gaussiana" :-). Também parece que seu limite está incorreto: variáveis ​​que são constantemente 1 satisfazem as condições desta pergunta. Para a=b e t=a21 a probabilidade é 1 , mas seu limite é assintoticamente 2exp(ca)0 como a cresce grande.
whuber
A probabilidade de todas as variáveis ​​serem 1 diminui exponencialmente. Acho que não entendi o seu comentário. Para e a I limite indicado é bastante trivialmente verdadeiro como a probabilidade da soma é maior do que éa=bt=a21t212(a1)eln(2)c(a1/a)
user1189053
11
Sinto muito por um erro meu. Eu pensei que tinha mencionado uniformemente acima. Então p = 1/2 e podemos tomar a e b maior do que qualquer constante (se necessário) para a desigualdade de espera
user1189053
2
A menos que meus olhos estejam me enganando, você está considerando uma soma de produtos, não um produto de somas. :-)
cardeal

Respostas:

7

A relação algébrica

S=i,jxiyj=ixijyj

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=i=1axi(a,1/2)0aY=j=1byj0b

Xa=1ai=1axi,

de onde

S=abXaXb=abZab.

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.aXaSab

O próximo passo é perceber que

Zab=XaXb=12((Xa+Xb2)2(XaXb2)2)=12(U2V2).

é 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.UVZabK0(|z|)/πabtPra,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 limiteabS±abSababmax(a,b)SabaK0a distribuição é mostrada em preto. (Como é simétrico em torno de , , basta observar a cauda negativa.)S0Pr(S>t)=Pr(S<t)

Figuras

À 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 .a500|S|<ab

Insets


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 .abab1030010160

s <- function(a, b) {
  # Returns the distribution of S as a vector indexed by its support.
  products <- factor(as.vector(outer(seq(-a, a, by=2), seq(-b, b, by=2))))
  probs <- as.vector(outer(dbinom(0:a, a, 1/2), dbinom(0:b, b, 1/2)))
  tapply(probs, products, sum)
}

par(mfrow=c(2,3))
b.vec <- c(51, 101, 149, 201, 299, 501)
cols <- terrain.colors(length(b.vec)+1)
for (a in c(50, 100, 150, 200, 300, 500)) {
  plot(c(-sqrt(a*max(b.vec)),0), c(10^(-300), 1), type="n", log="y", 
       xlab="S/sqrt(ab)", ylab="CDF", main=paste(a))
  curve(besselK(abs(x), 0)/pi, lwd=2, add=TRUE)
  for (j in 1:length(b.vec)) {
    b <- b.vec[j]
    x <- s(a,b)
    n <- as.numeric(names(x))
    k <- n <= 0
    y <- cumsum(x[k])
    lines(n[k]/sqrt(a*b), y, col=cols[j], lwd=2)
  }
}
whuber
fonte
11
Muito bem feito! Pode-se obter uma forma exata para o cdf do produto de 2 normais normais .. para a cauda negativa, é 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. a=5,b=7
wolfies
11
@wolfies Sim, eu também obtive essa expressão: integra a cauda de . Como a distribuição exata se afasta dela nas caudas extremas, não parecia valer a pena levar adiante a análise dessa integral. O próximo passo lógico é uma análise mais criteriosa das caudas, o que significa ir além da aproximação Normal. K0
whuber
3

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,...,abZXY

MZ(t)=E[ezt]=12et+12et=cosh(t)

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 éZW=Z1+Z2{2,0,2}{1/4,1/2,1/4}

MW(t)=E[e(z1+z2)t]=14e2t+12+14e2t==14(e2t+1)+14(e2t+1)=142etcosh(t)+142etcosh(t)=cosh(t)cosh(t)=MZ1(t)MZ2(t)

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 cadeia Zij=XiYj

P[Zab,...,Z11]=P[ZabZa,b1,...,Z11]...P[Z13Z12,Z11]P[Z12Z11]P[Z11]

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[Z12Z11]=P[Z12]
P[Z13,Z12Z11]Z13Z12Z11

P[Z13Z12,Z11]=P[Z13Z11]=P[Z13]

P[Z13Z12,Z11]P[Z12Z11]P[Z11]=P[Z13,Z12,Z11]=P[Z13]P[Z12]P[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

Sab=k=1abZk

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=evenS[ab,...,ab]ab=oddS[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
mZ1SS{ab2m;mZ+{0};mab}mSZ{Z1=z1,...,Zab=zab}S

P(S=ab2m)=(abm)12ab,0mab

Definindo , e número ímpar por construção, e o elemento típico do suporte de , temossab2mS

P(S=s)=(ababs2)12ab

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=oddS|S|

P(|S|=|s|)=(ababs2)12ab1

Então a função de distribuição é

P(|S||s|)=12ab11is,iodd(ababi2)

Portanto, para qualquer real , , obtemos a probabilidade requerida t1t<ab

P(|S|>t)=1P(|S|t)=112ab11it,iodd(ababi2)

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.5i9

Alecos Papadopoulos
fonte
O número de valores negativos em deve ser par . Portanto, essas quatro variáveis ​​aleatórias (presumo que sejam quatro dos seus - a notação não é clara) não são independentes. (X1Y1,X1Y2,X2Y1,X2Y2)Z
whuber
@whuber Obrigado. O problema (meu problema) é que continuo tendo independência em qualquer exemplo específico que eu resolvo. Vou trabalhar as quatro variáveis ​​específicas que você escreve.
Alecos Papadopoulos
Sim, é complicado porque distintos são independentes por pares e (acredito) quaisquer três distintos também são independentes. (I upvoted sua resposta por causa de seu ataque criativo sobre o problema e espero que eu esteja enganado na minha avaliação da falta de independência!)ZZ
whuber
@whuber Mais uma vez obrigado whuber, isso é realmente favorável. Estou pensando, o que precisamos para que a derivação da distribuição de seja válida é que todos os eventos sejam equiprobáveis. É possível uma propriedade assim, enquanto a independência conjunta falha? Quero dizer, a independência conjunta é suficiente para a equiprobabilidade, mas também é necessário? S{k=1abZk}
Alecos Papadopoulos
Receio não entender sua notação, que parece se referir a uma interseção de variáveis ​​aleatórias (o que isso possa significar).
whuber
3

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)

Sab=k=1abZkwhereZk=XiYj

"... 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:

  • pontos triângulo: simulação de Monte Carlo do pmf de dado eSa=5b=7
  • pontos redondos: simulação de Monte Carlo de uma caminhada aleatória simples com passosn=35

insira a descrição da imagem aqui

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 :Sab

 SumAB[a_, b_] :=  Outer[Times, RandomChoice[{-1, 1}, a], RandomChoice[{-1, 1}, b]] 
                         // Flatten // Total 

Em seguida, 500.000 tais caminhos, dizer quando e , pode ser gerada com:a=5b=7

 data57 = Table[SumAB[5, 7], {500000}];

O domínio de apoio para esta combinação de e é:ab

{-35, -25, -21, -15, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 15, 21, 25, 35}
wolfies
fonte
11
+1 Uma simulação (ou algum exemplo concreto) há muito tempo é necessária para nos fornecer uma referência para análises adicionais. Sua simulação pode ser muito mais eficiente (cerca de 25 vezes mais rápida) observando que como . Isso explica imediatamente por que nenhum valor primo suficientemente grande pode aparecer no seu gráfico de triângulos - e demonstra à força que não pode ter uma distribuição de "passeio aleatório" (Binomial em escala). S(ixi)(jyj)S
whuber
11
Em vez de simular, você pode obter rapidamente a resposta exata (para ae bmenos de 1000, pelo menos) como rademacher[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, digamos s[100,211],.
whuber
@ whuber re primeiro comentário - sua fatoração é super legal! :) No meu Mac, usando: ......... WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]]... é duas vezes mais rápido que a Outerabordagem. Curioso para saber qual código você está usando? [Ambas as abordagens podem, é claro, ser feita mais rapidamente utilizando ParallelTable, etc]
wolfies
Tente isto: 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 hora Tally[sum[500000, 5, 7]]. Para Raficianodos, 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").
whuber
@whuber - re comment2 - pmf exato: então você tem , onde cada soma da Rademacher é um binomial e, portanto, temos o produto de 2 binômios. Por que não escrever isso como resposta !? - é consideravelmente, puro, elegante e útil ...S=(iXi)(jYj)
wolfies