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ć . 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)?
źródło
Odpowiedzi:
Relacja algebraiczna
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=∑ai=1xi (a,1/2) 0 a Y=∑bj=1yj 0 b
skąd
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.a Xa S ab−−√
Następnym krokiem jest zauważenie tego
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.U V Zab K0(|z|)/π a b t Pra,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 ograniczeniaa b S ±ab S ab−−√ abmax(a,b)−−−−−−−−−−√ S a b a K0 dystrybucja jest pokazana na czarno. (Ponieważ jest symetryczny wokół , , więc wystarczy spojrzeć na ujemny ogon.)S 0 Pr(S>t)=Pr(−S<−t)
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 .a 500 |S|<ab√
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 .a b 10−300 10−160
a
b
źródło
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.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,...,ab Z X Y
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 jestZ W=Z1+Z2 {−2,0,2} {1/4,1/2,1/4}
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ńcuchaZij=XiYj
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[Z12∣Z11]=P[Z12]
P[Z13,Z12∣Z11] Z13 Z12 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
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=even S [−ab,...,ab] ab=odd S [−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
m Z −1 S S∈{ab−2m;m∈Z+∪{0};m≤ab} m S Z {Z1=z1,...,Zab=zab} S
Definiując i liczbę nieparzystą według konstrukcji oraz typowy element podparcia , mamys≡ab−2m S
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=odd S |S|
Zatem funkcja dystrybucji to
Dlatego dla dowolnego rzeczywistego , otrzymujemy wymagane prawdopodobieństwot 1≤t<ab
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.5 i 9
źródło
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)
„… 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:
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 :S a b
Następnie 500000 takie ścieżki, np gdy , a , można wytworzyć z:a=5 b=7
Dziedziną wsparcia dla tej kombinacji i jest:a b
źródło
a
ib
zarówno mniej niż 1000, w każdym razie) jakorademacher[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, powiedzmys[100,211]
.WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]]
... jest dwa razy szybszy niżOuter
podejście. Ciekawy, jakiego kodu używasz? [Oba podejścia można oczywiście przyspieszyć przy użyciuParallelTable
itp.]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 czasTally[sum[500000, 5, 7]]
. DlaR
aficianodos 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")
.