1
$\begingroup$

If I'm doing a simulation with $n$ trials, each with probability $p$, a quick way to select the successful trials is to choose a binomially distributed random number. Then randomly choose that many trials to be successful.

But what if the probabilities are distinct? I want to - as efficiently as possible - do a simulation in which I select the events that are successful. So I'd like to avoid generating a random number for each trial.

To make things specific, let's assume I've got $p_1$, ... , $p_{20}$ with each probability being some small number, say no bigger than $0.05$. I'd like to generate the successes without 20 coin flips.

Any suggestions? I know one way to do it if I order the probabilities and use a rejection sampling approach, but it would be great if I could avoid the cost of ordering them.

$\endgroup$
7
  • $\begingroup$ For the specific example: generate 20 random numbers at once and compare with the vector of success probabilities. As an example, consider the following Matlab implementation. Let p be the row vector of success probabilities. Then the row vector of booleans s = rand(1,20) > p gives you the events that are successful. $\endgroup$
    – Ritz
    Jan 16, 2015 at 12:34
  • $\begingroup$ I'm trying to avoid generating 20 random numbers. In this case, on average, less than one event is successful, so it seems likely that there may be something intelligent which would avoid generating so many numbers $\endgroup$
    – Joel
    Jan 16, 2015 at 13:18
  • $\begingroup$ Its not clear to me why your approach for homogenous probabilities is substantially faster than just simulating $n$ Bernoulli trials? Also, your criticism of Ritz's method could equally apply to homogenous probabilities. Rare event sampling is generally simulation intensive, unless you can come up with a clever importance sampling scheme. $\endgroup$
    – user76844
    Jan 16, 2015 at 15:18
  • $\begingroup$ The runtime for generating a number from a binomial distribution is roughly constant. The runtime of simulating n trials is O(n). $\endgroup$
    – Joel
    Jan 16, 2015 at 15:56
  • $\begingroup$ @Joel yes, but then what is the runtime for randomly selecting which trials are "successful" $\endgroup$
    – user76844
    Jan 16, 2015 at 16:25

2 Answers 2

Reset to default
0
$\begingroup$

Here's a suggestion:

Lets say you have $n$ trials ($T_i$), where $p_i = P(T_i=1), q_i=1-p_i$. Then there are $2^n$ combinations of successful trials (e.g., $(1,0,0,1)$ is one such combination if $n=4$).

  1. Set an $2^n$-vector, $\mathbf{P}:=(0,0,...0)$
  2. Let $I:=\{0, 1,2,...2^n-1\}$ be a set of integers
  3. For each $i\in I$, calculate its binary expansion as an n-vector $\mathbf{B}_i$ and
  4. Calculate $P(\mathbf{B}_i)=\prod\limits_{j=1}^n[p_iB_{ij}+q_i(1-B_{ij})]$
  5. Set $P_i=P(\mathbf{B}_i)$

After doing this, you will have a discrete probability distribution over $i\in I:\mathbf{P}(i)$

Now, to simulate successes, you just do the following:

  1. Let $X\sim \mathbf{P}$ and draw an integer from it.
  2. The successful trials will be the random vector $\mathbf{B}_X$

So, I've substituted some up-front computation to reduce your problem to a single-variable simulation.

$\endgroup$
0
$\begingroup$

You can efficiently generate the number of successful trials via the Poisson-Binomial distribution.

The PMF need only be calculated once for a given set of $P_i$, and variates generated by the method of your choice - using simple inverse sampling I generated >5 million / second on an old atom-based netbook.

Once you have the number of successes (or a sampling of them), you can use your alluded to method to select the actual members of the trial vector that were the successes. You might also have a peek at example 22.10.2 in Numerical Analysis for Statisticians by Lange.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .