# Variable probability Bernoulli outcomes – Fast and Slow

I am working on a project that requires the generation of Bernoulli outcomes. Typically, I would go about this using the built in sample() function like so:

```sample(1:0,n,prob=c(p,1-p),replace=TRUE)
```

This works great and is fast, even for large n. Problem is, I want to generate each sample with its own unique probability. Seems straight forward enough, I just wrapped the function and vectorized to allow the passing of a vector of p.

```binomial_sampler<-function(p){
return(sample(1:0,1,prob=c(p,1-p)))
}
bs<-Vectorize(binomial_sampler)
```

Naming this function bs() turned out to be rather prophetic. Nevertheless, I can call this function by passing my unique vector of outcome probabilities. And indeed I get the result I’m looking for.

```bs(my_p_vec)
```

Problem is, this turns out to be very slow. It would seem that there is quite a bit of overhead to calling sample() for one sample at a time. R’s RNGs are very fast for generating many iid samples, so I started thinking like my old c++ programming self and tried a different approach.

```Nbs<-function(p)
{
U<-runif(length(p),0,1)
outcomes<-U<p
return(outcomes)
}
```

I call the new version Nbs for “New Bernoulli Sampler”, or “Not Bull Shit”. And what a difference it made indeed!

```library(rbenchmark)
p<-runif(1000)
res <- benchmark(bs(p), Nbs(p))
print(res)
test replications elapsed relative user.self sys.self user.child sys.child
2 Nbs(p)          100   0.007        1     0.008    0.000          0         0
1  bs(p)          100   1.099      157     1.080    0.016          0         0
```

157x faster! Now that’s a speedup to write home about.

## 6 thoughts on “Variable probability Bernoulli outcomes – Fast and Slow”

1. edielivon says:

Just for kicks:

require(parallel)
require(doMC)
registerDoMC(2)
require(plyr)

bsply<-function(p){
laply(.data=p,
function(p)sample(1:0,1,prob=c(p,1-p)),
.parallel=T)}

test replications elapsed relative user.self sys.self user.child sys.child
2 Nbs(p) 100 0.006 1 0.005 0.000 0.000 0.000
1 bsply(p) 100 93.930 15655 66.293 5.631 7.975 12.794

2. Dason Kurkiewicz says:

This is a little bit slower than your proposed alternative but it seems natural to me to think of rbinom(length(p), 1, p) as another alternative.

• stevencarlislewalker says:

This is what I would have done too. This just makes Corey’s function even more interesting because it seems faster than what I would consider the standard base R thing (i.e. using rbinom).

3. dogg13 says:

The pedant in me can’t leave it alone…

Surely Nbs should be
{

return (outcomes)
}

Or maybe even
{

outcomes <- ifelse(U < p, 1, 0)
return(outcomes)
}

Cheers.

• Corey Chivers says:

Yes, you’re totally right – that was a typo. Did I mention that this software comes with ABSOLUTELY NO WARRANTY. Thanks for the catch 😉