7
$\begingroup$

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

Improve this question
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Daniel Lichtblau
    Mar 13, 2012 at 20:53
  • 2
    $\begingroup$ 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. $\endgroup$
    – rm -rf ♦
    Mar 13, 2012 at 21:12
  • $\begingroup$ 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? $\endgroup$
    – Sjoerd C. de Vries
    Mar 13, 2012 at 23:10

1 Answer 1

Reset to default
8
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Mar 13, 2012 at 21:25
  • $\begingroup$ 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? $\endgroup$
    – Sjoerd C. de Vries
    Mar 13, 2012 at 23:13
  • $\begingroup$ 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. $\endgroup$
    – whuber
    Mar 13, 2012 at 23:16
  • $\begingroup$ 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. $\endgroup$
    – Kevin O'Bryant
    Mar 15, 2012 at 1:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.