Summe der Produkte von Rademacher-Zufallsvariablen

9

Sei unabhängige Zufallsvariablen, die Werte oder mit einer Wahrscheinlichkeit von jeweils 0,5 . Betrachten Sie die Summe . Ich möchte die Wahrscheinlichkeit nach oben begrenzen . Die beste Grenze, die ich ist wobei c eine universelle Konstante ist. Dies wird erreicht, indem die Wahrscheinlichkeit Pr (| x_1 + \ Punkte + x_n | <\ sqrt {t}) und Pr (| y_1 + \ Punkte + y_n | <\ sqrt {t}) durch Anwendung einfacher Chernoff-Grenzen niedriger begrenzt werden. Kann ich hoffen, etwas zu bekommen, das wesentlich besser ist als diese Grenze? Für den Anfang kann ich wenigstens bekommenx1xa,y1yb+11S=i,jxi×yjP(|S|>t)2ectmax(a,b)cPr(|x1++xn|<t)Pr(|y1++yn|<t)ectab . Wenn ich subgaußsche Schwänze bekommen kann, wäre das wahrscheinlich das Beste, aber können wir das erwarten (ich denke nicht, aber ich kann mir kein Argument vorstellen)?

user1189053
quelle
Haben Sie darüber nachgedacht, den Chernoff direkt an S zu binden S? Möglicherweise können Sie etwas mit
E[exp(λS]=E[λijXiYj]=E[λ(iXi)(jYj)]
Dilip Sarwate
Es gibt eine offensichtliche Verbesserung Ihrer Grenze für t>ab , denn dann muss die Wahrscheinlichkeit Null sein. Es scheint mir, dass das ein "sub-Gaußscher" Schwanz ist :-). Es scheint auch, dass Ihre Bindung falsch ist: Variablen, die ständig 1 erfüllen die Bedingungen dieser Frage. Für a=b und t=a21 die Wahrscheinlichkeit 1 aber Ihre Grenze ist asymptotisch 2exp(ca)0 wenn a groß wird.
whuber
Die Wahrscheinlichkeit, dass alle Variablen 1 sind, sinkt exponentiell. Ich glaube nicht, dass ich Ihren Kommentar verstehe. Für und die von mir angegebene Grenze ganz trivial wahr, da die Wahrscheinlichkeit, dass die Summe größer als ist,a=bt=a21t212(a1)eln(2)c(a1/a)
user1189053
1
Es tut mir wirklich leid wegen eines Fehlers von mir. Ich dachte, ich hätte oben einheitlich erwähnt. Also ist p = 1/2 und wir können a und b größer als jede Konstante (falls erforderlich) nehmen, damit die Ungleichung gilt
user1189053
2
Wenn meine Augen mich nicht täuschen, betrachten Sie eine Summe von Produkten, kein Produkt von Summen. :-)
Kardinal

Antworten:

7

Die algebraische Beziehung

S=i,jxiyj=ixijyj

zeigt als Produkt zweier unabhängiger Summen. Da und unabhängige Bernoulli -Variaten sind, ist eine Binomialvariable , die wurde verdoppelt und verschoben. Daher ist sein Mittelwert und seine Varianz ist . In ähnlicher Weise hat einen Mittelwert von und eine Varianz von . Lassen Sie uns sie jetzt durch Definieren standardisierenS(xi+1)/2(yj+1)/2(1/2)X=i=1axi(a,1/2)0aY=j=1byj0b

Xa=1ai=1axi,

woher

S=abXaXb=abZab.

Um ein hohes (und quantifizierbares) Grad an Genauigkeit, wie wächst großer nähert sich die Standardnormalverteilung. Lassen Sie uns daher als mal das Produkt zweier Standardnormalen approximieren.aXaSab

Der nächste Schritt ist, das zu bemerken

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

ist ein Vielfaches der Differenz der Quadrate der unabhängigen Standard - Normalvariablen und . Die Verteilung von kann analytisch berechnet werden (durch Invertieren der charakteristischen Funktion ): Sein PDF ist proportional zur Bessel-Funktion der Ordnung Null, . Da diese Funktion exponentielle Schwänze hat, schließen wir sofort, dass es für großes und und festes keine bessere Annäherung an als in der Frage angegeben.UVZabK0(|z|)/πabtPra,b(S>t)

Es bleibt Raum für Verbesserungen, wenn eines (mindestens) von und nicht groß ist oder an Punkten im Schwanz von nahe . Direkte Berechnungen der Verteilung von zeigen eine gekrümmte Verjüngung der Schwanzwahrscheinlichkeiten an Punkten, die viel größer als , ungefähr jenseits von . Diese logarithmisch linearen Diagramme der CDF von für verschiedene Werte von (in den Titeln angegeben) und (ungefähr in den gleichen Werten wie , unterschieden durch Farbe in jedem Diagramm) zeigen, was los ist. Als Referenz dient der Graph der BegrenzungabS±abSababmax(a,b)SabaK0Verteilung wird in schwarz angezeigt. (Da um symmetrisch ist , ist , so dass es ausreicht, den negativen Schwanz zu betrachten.)S0Pr(S>t)=Pr(S<t)

Zahlen

Wenn größer wird, wächst die CDF näher an der Referenzlinie.b

Die Charakterisierung und Quantifizierung dieser Krümmung würde eine genauere Analyse der normalen Annäherung an Binomialvariablen erfordern.

Die Qualität der Bessel-Funktionsnäherung wird in diesen vergrößerten Bereichen (in der oberen rechten Ecke jedes Diagramms) klarer. Wir sind schon ziemlich weit draußen. Obwohl die logarithmische vertikale Skala erhebliche Unterschiede verbergen kann, eindeutig durch die Zeit erreicht die Annäherung ist gut für .a500|S|<ab

Einschübe


R-Code zur Berechnung der Verteilung vonS

Die Ausführung des folgenden Vorgangs dauert einige Sekunden. (Es werden mehrere Millionen Wahrscheinlichkeiten für 36 Kombinationen von und berechnet .) Lassen Sie auf langsameren Maschinen die größeren ein oder zwei Werte von und weg und erhöhen Sie die untere Plotgrenze von auf etwa .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
quelle
1
Sehr schön gemacht! Man kann eine genaue Form für das cdf des Produkts von 2 Standardnormalen erhalten. Für den negativen Schwanz ist es 1/2 (1 + y BesselK[0,-y] StruveL[-1, y] - y BesselK[1,-y] StruveL[0, y]). Es wäre interessant zu sehen, wie: (a) die Grenze des OP funktioniert und (b) Ihre normale Näherung für den oben betrachteten Fall funktioniert, dh abgeleitet unter Verwendung der exakten diskreten pmf-Lösung. a=5,b=7
Wolfies
1
@wolfies Ja, ich habe auch diesen Ausdruck erhalten: Er integriert den Schwanz von . Da die genaue Verteilung in den extremen Schwänzen davon abweicht, schien es nicht sinnvoll, die Analyse dieses Integrals weiter zu führen. Der logische nächste Schritt ist eine differenziertere Analyse der Schwänze, was bedeutet, über die normale Näherung hinauszugehen. K0
whuber
3

Kommentar: Ich habe den Titel bearbeitet, um besser zu reflektieren, welche Art von Wohnmobilen in der Frage berücksichtigt werden. Jeder kann sie erneut bearbeiten.

Motivation: Ich denke, es besteht keine Notwendigkeit, sich mit einer Obergrenze zufrieden zu geben, wenn wir die Verteilung vonableiten können . ( UPDATE : Wir können Whubers Kommentare und Antworten nicht sehen).|Sab|

Bezeichne . Es ist leicht zu überprüfen, ob die gleiche Verteilung wie und . Die Momenterzeugungsfunktion istZk=XiYj,k=1,...,abZXY

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

Darüber hinaus sind die paarweise unabhängig: Die Variable (Indizes können natürlich beliebig sein) unterstützt mit entsprechenden Wahrscheinlichkeiten . Seine Momenterzeugungsfunktion istZW=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)

Ich werde versuchen zu vermuten, dass die vollständige Unabhängigkeit wie folgt gilt (ist es für die weiseren offensichtlich?): Bezeichnen Sie für diesen Teil . Dann ist nach der Kettenregel Zij=XiYj

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

Durch paarweise Unabhängigkeit haben wir . Betrachten Sie . und sind unabhängig von also haben wir die zweite Gleichheit durch paarweise Unabhängigkeit. Aber das impliziert dasP[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]

Usw. (glaube ich). ( UPDATE : Ich denke falsch . Unabhängigkeit gilt wahrscheinlich für jedes Triplett, aber nicht für den gesamten Haufen. Was folgt, ist nur die Ableitung der Verteilung eines einfachen zufälligen Spaziergangs und keine korrekte Antwort auf die Frage - siehe Wolfies und Whubers Antworten).

Wenn die volle Unabhängigkeit tatsächlich zutrifft, haben wir die Aufgabe, die Verteilung einer Summe der dichotomen rvs

Sab=k=1abZk

das sieht aus wie ein einfacher zufälliger Spaziergang , allerdings ohne die klare Interpretation des letzteren als Sequenz.

Wenn die Unterstützung von die geraden ganzen Zahlen in einschließlich Null, während wenn die Unterstützung von die ungeraden ganzen Zahlen in ohne Null. ab=evenS[ab,...,ab]ab=oddS[ab,...,ab]

Wir behandeln den Fall von . Bezeichne als die Anzahl der , die den Wert . Dann wird die Unterstützung der geschrieben werden . Für jeden gegebenen , wir einen eindeutigen Wert für erhalten . Darüber hinaus sind aufgrund symmetrischer Wahrscheinlichkeiten und Unabhängigkeit (oder nur Austauschbarkeit?) Alle möglichen gemeinsamen Realisierungen der Variablen wahrscheinlich. Also zählen wir und stellen fest, dass die Wahrscheinlichkeitsmassenfunktion von ist,ab=odd
mZ1SS{ab2m;mZ+{0};mab}mSZ{Z1=z1,...,Zab=zab}S

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

Definieren und ungerade Anzahl von Konstruktion und das typische Element der Unterstützung von , haben wirsab2mS

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

Wechseln zu, da, wenn , die Verteilung von um Null symmetrisch ist, ohne die Wahrscheinlichkeitsmasse Null zuzuweisen, und somit die Verteilung vonwird durch "Falten" des Dichtediagramms um die vertikale Achse erhalten, wodurch die Wahrscheinlichkeiten für positive Werte im Wesentlichen verdoppelt werden.|S|ab=oddS|S|

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

Dann ist die Verteilungsfunktion

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

Daher erhalten wir für jedes reelle , , die erforderliche Wahrscheinlichkeit t1t<ab

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

Beachten Sie, dass die Angabe garantiert, dass die Summe nur bis zu den Werten läuft, die in der Unterstützung von- Wenn wir zum Beispiel , läuft immer noch bis , da es zwangsläufig ungerade ist, zusätzlich zu einer ganzen Zahl.i=odd|S|t=10.5i9

Alecos Papadopoulos
quelle
Die Anzahl der negativen Werte in muss gerade sein . Daher sind diese vier Zufallsvariablen (ich nehme an, sie sind vier Ihrer - die Notation ist unklar) nicht unabhängig. (X1Y1,X1Y2,X2Y1,X2Y2)Z
whuber
@whuber Danke. Das Problem (das heißt mein Problem) ist, dass ich in einem bestimmten Beispiel, das ich ausarbeite, immer wieder Unabhängigkeit erhalte. Ich werde die spezifischen vier Variablen bearbeiten, die Sie schreiben.
Alecos Papadopoulos
Ja, es ist schwierig, weil verschiedene paarweise unabhängig sind und (glaube ich) auch drei verschiedene unabhängig sind. (Ich habe Ihre Antwort wegen ihres kreativen Angriffs auf das Problem positiv bewertet und hoffe, dass ich mich in meiner Einschätzung des Mangels an Unabhängigkeit ZZ
irre
@whuber Nochmals vielen Dank whuber, das ist wirklich unterstützend. Ich denke, was wir brauchen, damit die Ableitung der Verteilung von gültig ist, ist, dass alle Ereignisse wahrscheinlich sind. Kann eine solche Immobilie gehalten werden, während die gemeinsame Unabhängigkeit versagt? Ich meine, die gemeinsame Unabhängigkeit reicht aus, um die Gleichwahrscheinlichkeit aufrechtzuerhalten, aber ist sie auch notwendig? S{k=1abZk}
Alecos Papadopoulos
Ich fürchte, ich verstehe Ihre Notation nicht, die sich auf einen Schnittpunkt von Zufallsvariablen zu beziehen scheint (was auch immer das bedeuten mag).
whuber
3

Keine Antwort, sondern ein Kommentar zu Alecos 'interessanter Antwort, der zu lang ist, um in ein Kommentarfeld zu passen.

Sei unabhängige Rademacher-Zufallsvariablen und sei unabhängige Rademacher-Zufallsvariable. Alecos bemerkt, dass:(X1,...,Xa)(Y1,...,Yb)

Sab=k=1abZkwhereZk=XiYj

"... sieht aus wie ein einfacher zufälliger Spaziergang ". Wenn es wie ein einfacher zufälliger Spaziergang wäre, wäre die Verteilung von symmetrisch "glockenförmig unimodal" um 0.S

Um zu veranschaulichen, dass es sich nicht um einen einfachen zufälligen Spaziergang handelt, finden Sie hier einen kurzen Monte-Carlo-Vergleich von:

  • Dreieckspunkte: Monte-Carlo-Simulation der pmf von bei undSa=5b=7
  • runde Punkte: Monte-Carlo-Simulation eines einfachen zufälligen Spaziergangs mit Schrittenn=35

Geben Sie hier die Bildbeschreibung ein

Offensichtlich ist kein einfacher zufälliger Spaziergang; Beachten Sie auch, dass S nicht auf alle geraden (oder ungeraden) ganzen Zahlen verteilt ist.S

Monte Carlo

Hier ist der Code (in Mathematica ), der verwendet wird, um eine einzelne Iteration der Summe zu erzeugen , wenn und :Sab

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

Dann können 500.000 solcher Pfade, beispielsweise wenn und , erzeugt werden mit:a=5b=7

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

Die Unterstützungsdomäne für diese Kombination von und ist:ab

{-35, -25, -21, -15, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 15, 21, 25, 35}
Wölfe
quelle
1
+1 Eine Simulation (oder ein solches konkretes Beispiel) ist seit langem erforderlich, um uns eine Referenz für die weitere Analyse zu geben. Ihre Simulation kann viel effizienter gestaltet werden (etwa 25-mal schneller), indem Sie feststellen, dass Faktoren . Dies erklärt sofort, warum in Ihrem Dreiecksdiagramm keine ausreichend großen Primwerte angezeigt werden können - und zeigt gewaltsam, dass keine "Random Walk" -Verteilung (skaliertes Binomial) haben kann. S(ixi)(jyj)S
whuber
1
Anstatt zu simulieren, können Sie schnell die genaue Antwort erhalten (für aund bbeide weniger als 1000 jedenfalls), wie 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]] z s[100,211].
whuber
@whuber re erster Kommentar - Ihre Faktorisierung ist super ordentlich! :) Auf meinem Mac mit: ......... WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]]... ist es doppelt so schnell wie der OuterAnsatz. Neugierig, welchen Code Sie verwenden? [Beide Ansätze können natürlich mit ParallelTableusw. schneller gemacht werden ]
Wolfies
Versuchen Sie Folgendes : 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]]. Dann Zeit Tally[sum[500000, 5, 7]]. Für RAficianodos macht Folgendes dasselbe und dauert nur 50% länger als 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 - genaue pmf: Sie haben also , wobei jede Summe von Rademachers ein Binomial ist, und so haben wir das Produkt von 2 Binomen. Warum nicht als Antwort aufschreiben? - Es ist hübsch, ordentlich, elegant und nützlich ...S=(iXi)(jYj)
Wolfies