# RandomVariate with a Discrete Distribution

Nature has provided me with a random variable $Z$ taking on the values $0, 1, 2, \ldots$, with probabilities $z_0, z_1, \cdots$. I can sample from the distribution of $Z$ reasonably efficiently (I have done so $2^{28}$ times), and so I have estimates $\hat{z_0},\hat{z_1},\ldots$ of $z_0,z_1,\ldots$.

I am interested in the number $f(z_0,z_1,\ldots,z_8)$, where $f$ is an explicit but involved function. Obviously, a good estimate is $f(\hat{z_0}, \hat{z_1}, \ldots,\hat{z_8})$, and a standard way to build a confidence interval is bootstrapping. To do this, I need to sample from the discrete distribution with probabilities $\hat{z_0}, \hat{z_1},\dots,\hat{z_8},q$, ($q$ chosen to make the probabilities add up to 1) and I need to sample $2^{28}$ times (I don't need the samples, just the number of times 1 comes up, 2 comes up, ...), and this process needs to be repeated, say, 1000 times.

That's a lot of sampling, and it's going way too slowly. The first law of fast Mathematica code is to use built-in functions. Is there a way to get RandomVariate to work with an arbitrary discrete distribution? Any other suggestions for rapidly sampling from an arbitrary distribution (again, I only need the Tally of the sampling, not the sample itself)?

-
Kevin, maybe instead change each $z_j$ into a smallish interval around the estimated value? Depending on the nature of $f$ you might get something not unreasonable. This assumes $f$ can handle Interval input, which may not be the case. –  Daniel Lichtblau Mar 13 '12 at 20:53
No time to answer right now, but perhaps EmpiricalDistribution is what you're looking for. It can be used with RandomVariate. If this works, I'll write an answer later in the day. –  rm -rf Mar 13 '12 at 21:12
You don't specify how you count the occurrences of each possible value of Z. Perhaps your bottleneck may be not the sampling but the counting? –  Sjoerd C. de Vries Mar 13 '12 at 23:10

This is a multinomial distribution. Obtain your bootstrap sample quickly as in this example:

z = RandomReal[{0, 1}, 10];
z = z / (Plus @@ z) (* Generate an example set of values for  z0^, ..., z8^, q *)
f = MultinomialDistribution[2^28, z];
Timing[RandomVariate[f, 1000];]


(1.342 seconds).

-
Assuming the observed counts of $0$ through $8$ are not too small (bigger than $30$ or so each)--you could do just fine using a Normal approximation to the distribution of the $\widehat{z}_i$ and assuming (slightly incorrectly) that they are independent. This would be about a thousand times faster. But that would be an important consideration only if you needed orders of magnitude more bootstrap draws. –  whuber Mar 13 '12 at 21:25
I assume that since the variance of a multinomial can be described analytically the confidence interval for f could be described analytically as well; so, no need for bootstrapping? –  Sjoerd C. de Vries Mar 13 '12 at 23:13
Good thought. It depends in part on the nature of $f$: we don't even know whether it's differentiable. Typically, one resorts to bootstrapping precisely when an analytic solution might not be trustworthy. –  whuber Mar 13 '12 at 23:16
Well, now I feel like a dope. The answer to my Mathematica question is "EmpiricalDistribution" (why is this not on the DiscreteDistribution page?), but this answer completely sidesteps what I thought the issue was. –  Kevin O'Bryant Mar 15 '12 at 1:31