# Distribution of product of sums of a set

Given a set $A$ of $n=kl$ positive real numbers, let $a_{ij}$ be a random number removed from the set $A$ (without replacement) for $i=1 \dots k$ and $j=1 \dots l$. What is the distribution and expected value of the following expression?

$$\prod_{i=1}^{k} \sum_{j=1}^{l} a_{ij}$$

Note: Using a set of 30 numbers for $l=3$ and $k=10$, I generated 10000 samples and computed the above formula for it. The following figure shows the distribution of the samples.

My question is that: how can I find the above distribution and expected value analytically?

Update

To makes the things clear, I past the code I used to generate the histogram here. This is a Mathematica code. The code generates a set of 30 random numbers between 0 and 1 (variable "nums"). Then it does the following 10000 times to generate the numbers for the smooth histogram: Shuffle the numbers in "nums" and assign it to "rands"; split the numbers in "rand" into 10 groups of 3 numbers; compute the sum of each 3 numbers; multiply the 10 summations to have the result.

SmoothHistogram@Module[{nums = RandomReal[1, {3*10}], rands},
Table[
rands = RandomSample@nums;
Times @@ Total /@ Table[rands[[{i, i + 1, i + 2}]], {i, 1, Length[nums], 3}],
{10000}
]
]

• Are the numbers in A randomly selected from the positive reals? May 30, 2012 at 10:04
• $A$ is a fixed given set. For instance $\{3,2,1,5,6,7\}$. But the distribution of the numbers is unknown. May 30, 2012 at 10:13
• The definition of the random variable is still unclear. Are you saying that the random variable is the product of the sum where one aij is missing from the sum and the choice of the missing number is randomly selected as the one left out? May 30, 2012 at 11:20
• @Michael I thought the question was clear. My interpretation of this random variable is operationally expressed in a line of R code in my reply (to simulate a set x of its values): one draws $l$ groups of $k$ values from a set $A$, without replacement, and forms the product of the sums of those values.
– whuber
May 30, 2012 at 16:13
• @whuber I disagree. It sounds to me like A is a fixed finite set of real numbers. You take one number away from the set and then compute the product of sums base on the remaining n-1=kl-1 aij. But your interpretation seems to jive more with the simulated example that the OP gave. Do you really think the question is clear? May 30, 2012 at 16:28

Provided $(a_{ij})$ is not "too bad," for large $k$ and large $l$ the distribution is approximately lognormal.

This follows because the sums will be approximately normal (by virtue of the Central Limit Theorem) and (for large $k$) almost iid. The product of a large number of iid (or almost iid) random variables is the exponential of the sum of their logs. Provided the log has a variance (which it must in this case, because all variables are bounded), the CLT again implies the sum of the logs is approximately normal, whence the exponential produces (by definition) a lognormal distribution.

(In truth the result will be a little more negatively skewed than a lognormal, because the logs of the sums will tend--in general--to be negatively skewed, whence so will the sum of those logs. In some cases, therefore, the histogram of simulated values may actually have a negative skew and not appear to be lognormal at all. This only means $k$ is not yet large enough for the CLT approximation for the sum of logs to be very good.)

By varying the elements of $(a_{ij})$, you can easily control two aspects of this (asymptotic) distribution: the shape and scale (corresponding to the spread and location of the logarithms). For instance, the histogram will look more like a classical (high coefficient of variation) lognormal distribution when $\left(a_{ij}\right)$ is highly skewed, as in this R simulation (again using $k=10$ and $l=3$):

a <- rgamma(3*10, shape=1, scale=1)
hist(replicate(10000, prod(apply(matrix(sample(a), ncol=10), 2, sum))))


Still assuming $(a_{ij})$ isn't "bad" and that $l$ and $k$ are large, the expectation of this product is approximately the expectation of the sums (because they are close to independent), which in turn are $l$ times the mean of $(a_{ij})$. Therefore, as a crude estimate, the expectation is approximately $l^k$ times the $l^\text{th}$ power of the mean of $(a_{ij})$. By virtue of the foregoing parenthetical skewness considerations, this would tend to be an underestimate.

For instance, I performed the preceding simulation for $16*16$ Gamma($0.1$) variates. The mean of $(a_{ij})$ was $0.117$; $1/1$6 of the $16^\text{th}$ root of the mean of x turned out to be $0.097$: reasonably close to the mean of $(a_{ij})$.

Finally, a "bad" set of values $(a_{ij})$ may have large clusters or extreme outliers: both of these throw the CLT approximations off. A "good" set of values will have a smoothly varying histogram and few or no outliers.

• Thank you for your nice answer. Some of the points is still unclear to me. My knowledge in probability theory is very basic, so some of the points might be obvious to others. 1) how can I estimate the parameters of the lognormal function? The mean is estimated by the expected value as you described, right? But how about stddev? 2) What can we say about smaller values of $l$ and $k$? Is there a general straightforward way to find the distribution in such cases? What if $l$ is small but $k$ is large? May 31, 2012 at 1:03
• (1) For some estimation methods see the Wikipedia article. (2) When either or both $l$ or $k$ are small, all bets are off: there will be no universal analytical expression. Deeper analysis of the cumulant generating function of the empirical distribution function of $(a_{ij})$ might produce something, but that would take a lot of work. For a practical application, why not just simulate the distribution? All it takes is the second line of R code in this reply.
– whuber
May 31, 2012 at 13:16

If my assumptions are right then contrary to the smooth curve shown above the distribution is actually discrete. The probability is P[X=Lrs] =1/(kl) for each r=1,2,...,k and s=1,2,...,l Lrs = ∏ ∑ aij summed over all i not equal to r and j not equal to s. The expectation is ∑∑Lij/(kl) for i=1,2,...k and j=1,2,...,l. The height of each distinct Lrs will be 1/(kl) and m/(kl) for values of Lrs that are common to m different combinations of r and s. To get a shape like the one shown above some values of m should be relatively large near the center of the distribution. The drop to 0 is more abrupt than what is shown for the distribution above. The probability mass function drops by multiples of 1/(kl).