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

 A: 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.
A: 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).
