Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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 1

up vote 7 down vote accepted

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

share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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