# Sum of Products of Rademacher random variables

Let $x_1 \ldots x_a,y_1 \ldots y_b$ be independent random variables taking values $+1$ or $-1$ with probability 0.5 each. Consider the sum $S = \sum_{i,j} x_i\times y_j$. I wish to upper bound the probability $P(|S| > t)$. The best bound I have right now is $2e^{-\frac{ct}{\max(a,b)}}$ where $c$ is a universal constant. This is achieved by lower bounding the probability $Pr(|x_1 + \dots + x_n|<\sqrt{t})$ and $Pr(|y_1 + \dots + y_n|<\sqrt{t})$ by application of simple Chernoff bounds. Can I hope to get something that is significantly better than this bound? For starters can I at least get $e^{-c\frac{t}{\sqrt{ab}}}$. If I can get sub-gaussian tails that would probably be the best but can we expect that (I don't think so but can't think of an argument)?

• Have you considered applying the Chernoff bound directly to $S$? You might be able to do something with $$E[\exp(\lambda S]=E\left[\lambda \sum_i\sum_j X_iY_j\right]=E\left[\lambda\left(\sum_i X_i\right)\left(\sum_j Y_j\right)\right]$$ Jul 9 '13 at 13:47
• There is an obvious improvement in your bound for $t \gt ab$, for then the probability must be zero. It seems to me that's a "sub-Gaussian" tail :-). It also seems your bound is incorrect: variables that are constantly $1$ satisfy the conditions of this question. For $a=b$ and $t=a^2-1$ the probability is $1$ but your bound is asymptotically $2\exp(-ca) \to 0$ as $a$ grows large.
– whuber
Jul 9 '13 at 13:48
• The probability of all the variables being 1 goes down exponentially. I dont think I understand your comment. For $a=b$ and $t=a^2 -1$ the bound I stated is quite trivially true as the probability of the sum being greater than $t^2 - 1$ is $2^{-(a-1)} \leq e^{-ln(2)c(a-1/a)}$ Jul 9 '13 at 15:56
• I am really sorry about a mistake of mine. I thought i had mentioned uniformally above. So p = 1/2 and we can take a and b larger than any constant (if needed) for the inequality to hold Jul 9 '13 at 20:32
• Unless my eyes are deceiving me, you are considering a sum of products, not a product of sums. :-) Jul 10 '13 at 1:30

The algebraic relation

$$S = \sum_{i,j} x_i y_j = \sum_i x_i \sum_j y_j$$

exhibits $S$ as the product of two independent sums. Because $(x_i+1)/2$ and $(y_j+1)/2$ are independent Bernoulli$(1/2)$ variates, $X=\sum_{i=1}^a x_i$ is a Binomial$(a, 1/2)$ variable which has been doubled and shifted. Therefore its mean is $0$ and its variance is $a$. Similarly $Y=\sum_{j=1}^b y_j$ has a mean of $0$ and variance of $b$. Let's standardize them right now by defining

$$X_a = \frac{1}{\sqrt a} \sum_{i=1}^a x_i,$$

whence

$$S = \sqrt{ab} X_a X_b = \sqrt{ab}Z_{ab}.$$

To a high (and quantifiable) degree of accuracy, as $a$ grows large $X_a$ approaches the standard Normal distribution. Let us therefore approximate $S$ as $\sqrt{ab}$ times the product of two standard normals.

The next step is to notice that

$$Z_{ab} = X_aX_b = \frac{1}{2}\left(\left(\frac{X_a+X_b}{\sqrt 2}\right)^2 - \left(\frac{X_a-X_b}{\sqrt 2}\right)^2 \right) = \frac{1}{2}\left(U^2 - V^2\right).$$

is a multiple of the difference of the squares of independent standard Normal variables $U$ and $V$. The distribution of $Z_{ab}$ can be computed analytically (by inverting the characteristic function): its pdf is proportional to the Bessel function of order zero, $K_0(|z|)/\pi$. Because this function has exponential tails, we conclude immediately that for large $a$ and $b$ and fixed $t$, there is no better approximation to ${\Pr}_{a,b}(S \gt t)$ than given in the question.

There remains some room for improvement when one (at least) of $a$ and $b$ is not large or at points in the tail of $S$ close to $\pm a b$. Direct calculations of the distribution of $S$ show a curved tapering off of the tail probabilities at points much larger than $\sqrt{ab}$, roughly beyond $\sqrt{ab\max(a,b)}$. These log-linear plots of the CDF of $S$ for various values of $a$ (given in the titles) and $b$ (ranging roughly over the same values as $a$, distinguished by color in each plot) show what's going on. For reference, the graph of the limiting $K_0$ distribution is shown in black. (Because $S$ is symmetric around $0$, $\Pr(S \gt t) = \Pr(-S \lt -t)$, so it suffices to look at the negative tail.)

As $b$ grows larger, the CDF grows closer to the reference line.

Characterizing and quantifying this curvature would require a finer analysis of the Normal approximation to Binomial variates.

The quality of the Bessel function approximation becomes clearer in these magnified portions (of the upper right corner of each plot). We're already pretty far out into the tails. Although the logarithmic vertical scale can hide substantial differences, clearly by the time $a$ has reached $500$ the approximation is good for $|S| \lt a\sqrt{b}$.

### R Code to Calculate the Distribution of $S$

The following will take a few seconds to execute. (It computes several million probabilities for 36 combinations of $a$ and $b$.) On slower machines, omit the larger one or two values of a and b and increase the lower plotting limit from $10^{-300}$ to around $10^{-160}$.

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))
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)
}
}

• Very nicely done! One can obtain an exact form for the cdf of the product of 2 standard Normals .. for the negative tail, it is 1/2 (1 + y BesselK[0,-y] StruveL[-1, y] - y BesselK[1,-y] StruveL[0, y]). It would be interesting to see how: (a) the OP's bound performs, and (b) your Normal approximation performs, for the case we were looking at above, i.e. $a=5, b= 7$ derived using the exact pmf discrete solution. Apr 17 '14 at 20:48
• @wolfies Yes, I obtained that expression too: it integrates the tail of $K_0$. Because the exact distribution does depart from it in the extreme tails, it did not seem worthwhile carrying the analysis of that integral any further. The logical next step is a more discerning analysis of the tails, which means going beyond the Normal approximation.
– whuber
Apr 18 '14 at 0:01

Comment: I edited the title in an attempt to reflect better what kind of r.v.'s are considered in the question. Anyone feel free to re-edit.

Motivation: I guess there is no need to settle for an upper bound, if we can derive the distribution of $|S_{ab}|$. (UPDATE: We can't -see Whuber's comments and answer).

Denote $Z_k = X_iY_j,\;\; k=1,...,ab$. It is easy to verify that $Z$'s have the same distribution as the $X$'s and the $Y$'s. The moment generating function is

$$M_Z(t) = E[e^{zt}]=\frac 12e^{-t}+\frac 12e^t = \cosh(t)$$

Moreover the $Z$'s are, to begin with, pair-wise independent: The variable $W = Z_1+Z_2$ (indices can be any of course), has support $\{-2,0,2\}$ with corresponding probabilities $\{1/4,1/2,1/4\}$. Its moment generating function is

$$M_{W}(t) = E[e^{(z_1+z_2)t}] = \frac 14e^{-2t}+\frac 12 +\frac 14e^{2t}=\\ =\frac 14(e^{-2t}+1)+ \frac 14(e^{2t}+1) = \frac 14 2e^{-t}\cosh(t)+\frac 14 2e^{t}\cosh(t)\\ =\cosh(t)\cdot \cosh(t) = M_{Z_1}(t)M_{Z_2}(t)$$

I will attempt to suspect that full independence holds, as follows (is it obvious to the wiser ones?): For this part, denote $Z_{ij}=X_iY_j$. Then by the chain rule $$P[Z_{ab},...,Z_{11}] = P[Z_{ab}\mid Z_{a,b-1},...,Z_{11}]\cdot ...\cdot P[Z_{13}\mid Z_{12},Z_{11}]\cdot P[Z_{12}\mid Z_{11}]\cdot P[Z_{11}]$$

By pair-wise independence we have $P[Z_{12}\mid Z_{11}] = P[Z_{12}]$.
Consider $P[Z_{13},Z_{12}\mid Z_{11}]$. $Z_{13}$ and $Z_{12}$ are independent conditional on $Z_{11}$ so we have $$P[Z_{13}\mid Z_{12},Z_{11}] = P[Z_{13}\mid Z_{11}] = P[Z_{13}]$$ the second equality by pair-wise independence. But this implies that

$$P[Z_{13}\mid Z_{12},Z_{11}]\cdot P[Z_{12}\mid Z_{11}]\cdot P[Z_{11}] = P[Z_{13},\,Z_{12},\,Z_{11}] = P[Z_{13}]\cdot P[Z_{12}]\cdot P[Z_{11}]$$

Etc (I think). ( UPDATE : I think wrong. Independence probably holds for any triplet, but not for the whole bunch. So what follows is just the derivation of the distribution of a simple random walk, and not a correct answer to the question - see Wolfies' and Whuber's answers).

If full independence does indeed hold, we have the task of deriving the distribution of a sum of i.i.d dichotomous r.v.'s $$S_{ab}=\sum_{k=1}^{ab}Z_k$$

which looks like a simple random walk, although without the clear interpretation of the latter as a sequence.

If $ab=even$ the support of $S$ will be the even integers in $[-ab,...,ab]$ including zero, while if $ab=odd$ the support of $S$ will be the odd integers in $[-ab,...,ab]$, without zero.

We treat the case of $ab=odd$.
Denote $m$ to be the number of $Z$'s taking the value $-1$. Then the support of $S$ can be written $S\in \{ab-2m;m\in \mathbb Z_+\cup\{0\};m\le ab\}$. For any given $m$, we obtain a unique value for $S$. Moreover, due to symmetric probabilities and independence (or just exchangeability?), all possible joint realizations of the $Z$-variables $\{Z_1=z_1,..., Z_{ab}=z_{ab}\}$ are equiprobable. So we count and we find that the probability mass function of $S$ is,

$$P(S=ab-2m)={ab \choose m}\cdot \frac 1{2^{ab}}, \qquad 0\le m\le ab$$

Defining $s\equiv ab-2m$, and odd number by construction, and the typical element of the support of $S$, we have

$$P(S=s)={ab \choose \frac{ab-s}{2}}\cdot \frac 1{2^{ab}}$$

Moving to $|S|$, since if $ab=odd$, the distribution of $S$ is symmetric around zero without allocating probability mass to zero, and so the distribution of $|S|$ is obtained by "folding" the density graph around the vertical axis, essentially doubling the probabilities for positive values,

$$P(|S|=|s|)={ab \choose \frac{ab-s}{2}}\cdot \frac 1{2^{ab-1}}$$

Then the distribution function is

$$P(|S|\le|s|)=\frac 1{2^{ab-1}}\sum_{1\le i\le s,\, i\,odd}{ab \choose \frac{ab-i}{2}}$$

Therefore, for any real $t$, $1\le t<ab$, we obtain the required probability $$P(|S|> t) = 1- P(|S|\le t) = 1-\frac 1{2^{ab-1}}\sum_{1\le i\le t,\, i\,odd}{ab \choose {\frac{ab-i}{2}}}$$

Note that the indication $i=odd$ guarantees that the sum will run only up to values included in the support of $|S|$ - for example, if we set $t=10.5$, still $i$ will run up to $9$, since it is constrained to be odd, on top of being an integer.

• The number of negative values in $(X_1Y_1,X_1Y_2,X_2Y_1,X_2Y_2)$ must be even. Therefore these four random variables (I presume they are four of your $Z$s--the notation is unclear) are not independent.
– whuber
Apr 16 '14 at 15:38
• @whuber Thanks. The problem (my problem, that is), is that I keep getting independence in any specific example I work out. I will work the specific four variables you write. Apr 16 '14 at 16:20
• Yes, it's tricky because distinct $Z$s are pairwise independent and (I believe) any three distinct $Z$s are independent, too. (I upvoted your answer because of its creative attack on the problem and I hope I am mistaken in my assessment of the lack of independence!)
– whuber
Apr 16 '14 at 17:49
• @whuber Thanks again whuber, that's really supportive. I am thinking, what we need in order for the derivation of the distribution of $S$ to be valid, is that all events $\{\cap_{k=1}^{ab}Z_k\}$ are equiprobable. Is it possible for such a property to hold, while joint independence fails? I mean, joint independence is sufficient for equiprobability to hold, but is it also necessary? Apr 16 '14 at 17:58
• I'm afraid I don't understand your notation, which appears to refer to an intersection of random variables (whatever that might mean).
– whuber
Apr 16 '14 at 18:00

Not an answer, but a comment on Alecos’s interesting answer that is too long to fit into a comment box.

Let $(X_1, ..., X_a)$ be independent Rademacher random variables, and let $(Y_1, ..., Y_b)$ be independent Rademacher random variables. Alecos notes that:

$$S_{ab}=\sum_{k=1}^{ab}Z_k \qquad \text{where} \qquad Z_k = X_i Y_j$$

"… looks like a simple random walk”. If it were like a simple random walk, then the distribution of $S$ would be symmetric 'bell-shaped unimodal' around 0.

To illustrate that it is not a simple random walk, here is a quick Monte Carlo comparison of:

• triangle dots: Monte Carlo simulation of the pmf of $S$ given $a = 5$ and $b = 7$
• round dots: Monte Carlo simulation of a simple random walk with $n = 35$ steps

Clearly, $S$ is not a simple random walk; also note that S is not distributed on all the even (or odd) integers.

Monte Carlo

Here is the code (in Mathematica) used to generate a single iteration of the sum $S$, given $a$ and $b$:

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


Then, 500,000 such paths, say when $a = 5$ and $b = 7$, can be generated with:

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


The domain of support for this combination of $a$ and $b$ is:

{-35, -25, -21, -15, -9, -7, -5, -3, -1, 1, 3, 5, 7, 9, 15, 21, 25, 35}

• +1 A simulation (or some such concrete example) has long been needed to give us a reference for further analysis. Your simulation can be made much more efficient (about 25 times faster) by noting that $S$ factors as $\left(\sum_i x_i\right)\left(\sum_j y_j\right)$. That immediately explains why no sufficiently large prime values can show up in your triangle chart--and forcibly demonstrates that $S$ cannot have a "random walk" (scaled Binomial) distribution.
– whuber
Apr 17 '14 at 15:47
• Instead of simulating you can quickly obtain the exact answer (for a and b both less than 1000, anyway) as 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]] Try it with, say, s[100,211].
– whuber
Apr 17 '14 at 16:20
• @whuber re first comment - your factorisation is super neat! :) On my Mac, using: ......... WHuberSumAB[a_, b_] := Total[RandomChoice[{-1, 1}, a]] * Total[RandomChoice[{-1, 1}, b]] ... it is twice as fast as the Outer approach. Curious as to what code you are using? [Both approaches can, of course, be made faster using ParallelTable, etc] Apr 17 '14 at 16:45
• Try this: 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]]. Then time Tally[sum[500000, 5, 7]]. For R aficianodos, the following does the same thing and takes only 50% longer than 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
Apr 17 '14 at 17:03
• @whuber - re comment2 - exact pmf: so you have $S = (\sum_{i}^{}X_i)(\sum_{j}^{}Y_j)$, where each sum of Rademacher's is a Binomial, and so we have the product of 2 Binomials. Why not write this up as an answer!? - it is pretty, neat, elegant and useful ... Apr 17 '14 at 19:34