# 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

2
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

## 1 Answer

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
lang-mma
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 21, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0}, "config": {"markdown_headings": true, "markdown_code": false, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "github"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "english_confidence": 0.9216378331184387, "perplexity_flag": "middle"}