Suma produktów zmiennych losowych Rademacher

9

Niech będą niezależnymi zmiennymi losowymi przyjmującymi wartości lub -1 z prawdopodobieństwem 0,5 każda. Rozważ sumę S = \ sum_ {i, j} x_i \ razy y_j . Chciałbym przekroczyć górną granicę prawdopodobieństwa P (| S |> t) . Najlepsza granica, jaką mam teraz, to 2e ^ {- \ frac {ct} {\ max (a, b)}}, gdzie c jest stałą uniwersalną. Osiąga się to poprzez niższe ograniczenie prawdopodobieństwa Pr (| x_1 + \ dots + x_n | <\ sqrt {t}) i Pr (| y_1 + \ dots + y_n | <\ sqrt {t}) poprzez zastosowanie prostych granic Chernoffa. Czy mogę mieć nadzieję, że dostanę coś znacznie lepszego niż to ograniczenie? Na początek mogę przynajmniej dostaćx1xa,y1yb+11S=i,jxi×yjP(|S|>t)2ectmax(a,b)cPr(|x1++xn|<t)Pr(|y1++yn|<t)ectab . Jeśli uda mi się uzyskać ogony sub-gaussowskie, które prawdopodobnie byłyby najlepsze, ale czy możemy się tego spodziewać (nie sądzę, ale nie mogę wymyślić kłótni)?

użytkownik1189053
źródło
Czy rozważałeś zastosowanie Chernoffa związanego bezpośrednio z S ? Możesz zrobić coś za pomocą
E[exp(λS]=E[λijXiYj]=E[λ(iXi)(jYj)]
Dilip Sarwate
Istnieje wyraźna poprawa w twoim limicie dla t>ab , ponieważ wtedy prawdopodobieństwo musi wynosić zero. Wydaje mi się, że to „sub-gaussowski” ogon :-). Wydaje się również, że twoja granica jest niepoprawna: zmienne, które stale 1 spełniają warunki tego pytania. W przypadku a=b oraz t=a21 jest prawdopodobieństwo, 1 , lecz związany jest swoją asymptotycznie 2exp(ca)0 w rośnie duża. a
whuber
Prawdopodobieństwo, że wszystkie zmienne będą równe 1, spada wykładniczo. Nie sądzę, że rozumiem twój komentarz. Dla i granica, którą podałem, jest dość trywialnie prawdziwa, ponieważ prawdopodobieństwo, że suma jest większa niż wynosia=bt=a21t212(a1)eln(2)c(a1/a)
user1189053
1
Naprawdę mi przykro z powodu mojego błędu. Myślałem, że wspominałem jednolicie powyżej. Więc p = 1/2 i możemy przyjąć aib większą niż jakakolwiek stała (w razie potrzeby) dla utrzymania nierówności
użytkownik1189053
2
O ile moje oczy mnie nie oszukują, rozważasz sumę produktów, a nie iloczyn sum. :-)
kardynał

Odpowiedzi:

7

Relacja algebraiczna

S=i,jxiyj=ixijyj

wykazuje jako iloczyn dwóch niezależnych sum. Ponieważ i są niezależnymi zmiennymi Bernoulliego , jest zmienną dwumianową , która został podwojony i przesunięty. Dlatego jego średnia wynosi , a jego wariancja jest . Podobnie ma średnią i wariancję . Standaryzujmy je teraz, definiującS(xi+1)/2(yj+1)/2(1/2)X=i=1axi(a,1/2)0aY=j=1byj0b

Xa=1ai=1axi,

skąd

S=abXaXb=abZab.

Do wysokiego (i policzalne) stopień dokładności, jak rośnie duża zbliża się do rozkładu normalnego. Przybliżmy zatem jako razy iloczyn dwóch standardowych normalnych.aXaSab

Następnym krokiem jest zauważenie tego

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

jest wielokrotnością różnicy kwadratów niezależnych zmiennych standardowych Zwykłe i . Rozkład można obliczyć analitycznie ( odwracając funkcję charakterystyczną ): jego pdf jest proporcjonalny do funkcji Bessela rzędu zero, . Ponieważ ta funkcja jest wykładniczy ogony natychmiast stwierdzić, że w przypadku dużych i i stałej , nie ma lepsze przybliżenie niż w pytaniu.UVZabK0(|z|)/πabtPra,b(S>t)

Pozostaje potrzeba tak, gdy jedna (co najmniej) z i jest niewielka lub w punktach ogona blisko . Bezpośrednie obliczenia rozkładu pokazują zakrzywione zwężenie prawdopodobieństwa ogona w punktach znacznie większych niż , mniej więcej powyżej . Te logarytmiczno-liniowe wykresy CDF dla dla różnych wartości (podanych w tytułach) (w przybliżeniu powyżej tych samych wartości co , rozróżnianych kolorem na każdym wykresie) pokazują, co się dzieje. Dla porównania wykres ograniczeniaabS±abSababmax(a,b)SabaK0dystrybucja jest pokazana na czarno. (Ponieważ jest symetryczny wokół , , więc wystarczy spojrzeć na ujemny ogon.)S0Pr(S>t)=Pr(S<t)

Ryciny

Gdy rośnie, CDF zbliża się do linii odniesienia.b

Charakterystyka i kwantyfikacja tej krzywizny wymagałaby dokładniejszej analizy przybliżenia normalnego do zmiennych dwumianowych.

Jakość przybliżenia funkcji Bessela staje się wyraźniejsza w tych powiększonych częściach (w prawym górnym rogu każdego wykresu). Jesteśmy już dość daleko w tyle. Chociaż logarytmiczna skala pionowa może ukryć znaczne różnice, wyraźnie do czasu, gdy osiągnie przybliżenie jest dobre dla .a500|S|<ab

Wypustki


Kod R do obliczenia rozkładuS

Wykonanie poniższych czynności potrwa kilka sekund. (Oblicza kilka milionów prawdopodobieństwa dla 36 kombinacji i ). Na dłuższych maszyn pominąć większy lub dwie wartości i i wzrost dolnej granicy kreślenia od około .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
źródło
1
Bardzo ładnie wykonane! Dokładną formę można uzyskać dla cdf produktu 2 standardowych normalnych .. dla negatywnego ogona, to jest 1/2 (1 + y BesselK[0,-y] StruveL[-1, y] - y BesselK[1,-y] StruveL[0, y]). Interesujące byłoby zobaczyć, jak: (a) wykonuje się operacja związana z OP, oraz (b) wykonuje się normalne przybliżenie, w przypadku, o którym mówiliśmy powyżej, tj. wyprowadzone przy użyciu dokładnego rozwiązania dyskretnego pmf. a=5,b=7
wilki
1
@wolfies Tak, też otrzymałem to wyrażenie: integruje ogon . Ponieważ dokładny rozkład odbiega od niego w skrajnych ogonach, nie warto wydawać dalszej analizy tej całki. Logicznym następnym krokiem jest bardziej szczegółowa analiza ogonów, co oznacza wyjście poza normalne przybliżenie. K0
whuber
3

Komentarz: Zredagowałem tytuł, aby lepiej odzwierciedlić, jakie rv są rozważane w pytaniu. Każdy może ponownie edytować.

Motywacja: Myślę, że nie ma potrzeby zadowolenia się górną granicą, jeśli możemy wyprowadzić rozkład. ( AKTUALIZACJA : Nie możemy zobaczyć komentarzy i odpowiedzi Whubera).|Sab|

Oznaczają . Jest to łatwe do sprawdzenia, że „y mają ten sam rozkład co ” S i „S. Funkcja generowania momentu toZk=XiYj,k=1,...,abZXY

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

Co więcej, są na początek niezależne parami: Zmienna (wskaźniki mogą być oczywiście dowolne), ma wsparcie z odpowiednimi prawdopodobieństwami . Jego funkcją generowania momentu jestZW=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)

Spróbuję podejrzewać, że zachowuje się pełna niezależność w następujący sposób (czy jest to oczywiste dla mądrzejszych?): W tej części . Następnie według reguły łańcucha Zij=XiYj

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

Dzięki niezależności parami mamy . Rozważmy . i są niezależne od więc mamy druga równość przez niezależność par. Ale to implikujeP[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]

Itp (myślę). ( AKTUALIZACJA : Wydaje mi się, że jest źle . Niepodległość prawdopodobnie dotyczy każdej trypletu, ale nie całej grupy. A więc po prostu wyprowadzenie rozkładu zwykłego losowego marszu, a nie poprawna odpowiedź na pytanie - patrz Wolfies i Odpowiedzi Whubera).

Jeśli rzeczywiście zachodzi pełna niezależność, mamy za zadanie wyprowadzić rozkład sumy Sid dychotomicznego rv iid d_hotomous rv

Sab=k=1abZk

który wygląda jak zwykły losowy spacer , choć bez jasnej interpretacji tego ostatniego jako sekwencji.

Jeśli wsparcie będzie parzystymi liczbami całkowitymi w włączając zero, podczas gdy jeśli wsparcie będzie nieparzystymi liczbami całkowitymi w , bez zera. ab=evenS[ab,...,ab]ab=oddS[ab,...,ab]

Traktujemy przypadek . Oznaczmy jako liczbę przyjmującą wartość . Następnie można napisać wsparcie dla . Dla danego , otrzymujemy unikalną wartość . Ponadto, ze względu na symetryczne prawdopodobieństw i niezależności (lub po prostu zamienności?), Wszystkich możliwych wspólnych realizacje -zmienne są jednakowo prawdopodobne. Liczymy więc i stwierdzamy, że funkcja masy prawdopodobieństwa jest,ab=odd
mZ1SS{ab2m;mZ+{0};mab}mSZ{Z1=z1,...,Zab=zab}S

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

Definiując i liczbę nieparzystą według konstrukcji oraz typowy element podparcia , mamysab2mS

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

Przeprowadzka do, ponieważ jeśli , rozkład jest symetryczny wokół zera bez przydzielania masy prawdopodobieństwa do zera, a zatem rozkładuzyskuje się przez „złożenie” wykresu gęstości wokół osi pionowej, zasadniczo podwajając prawdopodobieństwo dla wartości dodatnich,|S|ab=oddS|S|

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

Zatem funkcja dystrybucji to

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

Dlatego dla dowolnego rzeczywistego , otrzymujemy wymagane prawdopodobieństwo t1t<ab

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

Zauważ, że wskazanie gwarantuje, że suma będzie działać tylko do wartości zawartych w obsłudze- Na przykład, jeśli zestaw , wciąż będą działać do , ponieważ jest ona ograniczona być dziwne, na górze jest liczbą całkowitą.i=odd|S|t=10.5i9

Alecos Papadopoulos
źródło
Liczba wartości ujemnych w musi być parzysta . Dlatego te cztery zmienne losowe (zakładam, że są to cztery z Twoich - notacja jest niejasna) nie są niezależne. (X1Y1,X1Y2,X2Y1,X2Y2)Z
whuber
@whuber Thanks. Problem (to znaczy mój problem) polega na tym, że uzyskuję niezależność w każdym konkretnym przykładzie, który wypracowuję. Będę pracował nad konkretnymi czterema zmiennymi, które napiszesz.
Alecos Papadopoulos
Tak, jest to trudne, ponieważ różne są niezależne parami i (sądzę) dowolne trzy różne są również niezależne. (Poparłem twoją odpowiedź z powodu jej twórczego ataku na problem i mam nadzieję, że się mylę w mojej ocenie braku niezależności!)ZZ
whuber
@ whuber Jeszcze raz dziękuję, to naprawdę wspierające. Myślę, że potrzebujemy, aby wyprowadzenie rozkładu było ważne, aby wszystkie zdarzenia były możliwe. Czy taka własność jest w stanie utrzymać, a wspólna niezależność zawodzi? Chodzi mi o to, że wspólna niezależność jest wystarczająca, aby utrzymać równowagę, ale czy jest również konieczna? S{k=1abZk}
Alecos Papadopoulos
Obawiam się, że nie rozumiem waszej notacji, która wydaje się odnosić do przecięcia zmiennych losowych (cokolwiek to może znaczyć).
whuber
3

Nie odpowiedź, ale komentarz do interesującej odpowiedzi Alecos, która jest zbyt długa, aby zmieściła się w polu komentarza.

Niech będą niezależnymi zmiennymi losowymi Rademacher i niech będą niezależnymi zmiennymi losowymi Rademacher. Alecos zauważa, że:(X1,...,Xa)(Y1,...,Yb)

Sab=k=1abZkwhereZk=XiYj

„… Wygląda jak zwykły losowy spacer ”. Gdyby to był zwykły losowy spacer, wówczas rozkład byłby symetryczny „unimodalny” w kształcie dzwonu wokół 0.S

Aby zilustrować, że to nie proste błądzenia losowego, tutaj jest szybkie porównanie Monte Carlo:

  • Kropki trójkąt: symulacji Monte Carlo PMF z danym , aSa=5b=7
  • okrągłe kropki: symulacja Monte Carlo prostego losowego marszu z krokówn=35

wprowadź opis zdjęcia tutaj

Oczywiście nie jest zwykłym przypadkowym spacerem; zauważ także, że S nie jest rozłożone na wszystkie parzyste (lub nieparzyste) liczby całkowite.S

Monte Carlo

Oto kod (w Mathematica ) użyty do wygenerowania pojedynczej iteracji sumy , biorąc uwagę i :Sab

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

Następnie 500000 takie ścieżki, np gdy , a , można wytworzyć z:a=5b=7

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

Dziedziną wsparcia dla tej kombinacji i jest:ab

{-35, -25, -21, -15, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 15, 21, 25, 35}
wilki
źródło
1
+1 Od dawna potrzebna była symulacja (lub jakiś taki konkretny przykład), aby dać nam odniesienie do dalszej analizy. Twoja symulacja może być znacznie wydajniejsza (około 25 razy szybsza), zwracając uwagę, że czynniki . To natychmiast tłumaczy, dlaczego na wykresie trójkątów nie mogą pojawić się wystarczająco duże liczby pierwsze - i na siłę pokazuje, że nie może mieć rozkładu „losowego przejścia” (skalowanego dwumianowego). S(ixi)(jyj)S
whuber
1
Zamiast symulacji można szybko uzyskać dokładną odpowiedź (dla ai bzarówno mniej niż 1000, w każdym razie) jako 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]] Spróbuj z, powiedzmy s[100,211].
whuber
@ whuber ponownie pierwszy komentarz - Twoja faktoryzacja jest bardzo fajna! :) Na moim Macu, używając: ......... WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]]... jest dwa razy szybszy niż Outerpodejście. Ciekawy, jakiego kodu używasz? [Oba podejścia można oczywiście przyspieszyć przy użyciu ParallelTableitp.]
wilki
Spróbuj tego: 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]]. Potem czas Tally[sum[500000, 5, 7]]. Dla Raficianodos dodaje robi to samo i trwa tylko 50% dłużej niż 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 komentarz 2 - dokładna pmf: więc masz , gdzie każda suma Rademachera jest dwumianowa, a więc mamy iloczyn 2 dwumianów. Dlaczego nie napisać tego jako odpowiedzi !? - jest ładny, schludny, elegancki i użyteczny ...S=(iXi)(jYj)
wilki