# How to sample vector of independent Bernoulli variables efficiently?
**URL:** https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563
**Category:** Statistics
**Created:** [January 12, 2019, 5:30pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563 "2019-01-12T17:30:23Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Gregstrq](https://discourse.julialang.org/u/Gregstrq)
#### Post date: [January 12, 2019, 5:30pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563/1 "2019-01-12T17:30:24Z")
Hi, everyone.
I am trying to perform disorder averaging over different realizations of a lattice with some of the nodes randomly thrown away.
Specifically, consider the following model: we have a lattice (graph) with N nodes. Then each of the nodes is thrown away with probability p. Realization of the lattice is represented as a vector of N zeros and ones, where each of the components is a random variable drawn from Bernoulli distribution (probablity p for zero and (1-p) for one). Then I need to average some quantity over this random vectors.
I’ve tried to generate random vectors with a simple bruteforce method `vec = [rand()
### Author:  [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)
#### Post date: [January 12, 2019, 7:12pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563/2 "2019-01-12T19:12:38Z")
Have you tried `Distributions.jl`? You can use `Binomial(n,p)` with `n = 1` like this:
```julia
using Distributions
p = 0.5
N = 100
d = Binomial(1,p)
v = rand(d, N)
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 12, 2019, 7:29pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563/3 "2019-01-12T19:29:25Z")
> [@Gregstrq](#):
>
> I am trying to perform disorder averaging over different realizations of a lattice with some of the nodes randomly thrown away.
You might also want to look at the `randsubseq` function in the `Random` standard library, which allows you to efficiently find a random (Bernouilli) subsequence of a given array. It is _far_ more efficient than calling `rand() < p` for every element of the array if `p` is much smaller than 1. (Or conversely you can apply it to the inverse draw if `1-p` if small.)
You can find the implementation [here](https://github.com/JuliaLang/julia/blob/19fedfbb686502a9fbd29cbe086e1b5e9d945965/stdlib/Random/src/misc.jl#L83-L121). In particular, if you are only computing some _function_ of the subsequence, you might want to adapt the implementation of `randsubseq` to compute your function directly, without actually allocating/constructing the subsequence array.
> [@Gregstrq](#):
>
> I’ve tried to generate random vectors with a simple bruteforce method `vec = [rand() […] The computation for each of the vectors takes about 100 seconds,
How can the computation of `vec` take 100 seconds? Even for `N=10^9`, on my laptop this computation takes only 3 seconds. Remember [not to benchmark in global scope](https://docs.julialang.org/en/latest/manual/performance-tips/#Avoid-global-variables-1) — put all of your performance-critical code into functions.
---
### Author:  [@Gregstrq](https://discourse.julialang.org/u/Gregstrq)
#### Post date: [January 12, 2019, 10:14pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563/4 "2019-01-12T22:14:21Z")
I am sorry I wrote my post in a confusing way.
For each of the realisations of the lattice I solve dynamical equations of motion and extract correlation functions. Then I need to average this correlation functions over all realizations. Integration of the system of equations for a particular realization takes about 100 seconds. This time is much larger then the time I spend generating realization of a lattice and setting up a system of equations of motion. So, I don’t care how much time I spend generating a vector. The problem is how many realizations I need to average over so that the result is very close to the one obtained by averaging over the infinite number of realizations. So by efficiency I mean the speed with which the average (of correlation functions) over a finite set of random vectors converges (when we make the size of the set larger and larger) to the value obtained for infinite set. Is it possible to speed this convergence by a suitable choice of sampling algorithm? (may be I need to use some analogue of Metropolis algorithm for discrete distribution).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 12, 2019, 11:10pm UTC](https://discourse.julialang.org/t/how-to-sample-vector-of-independent-bernoulli-variables-efficiently/19563/5 "2019-01-12T23:10:06Z")
> [@Gregstrq](#):
>
> So by efficiency I mean the speed with which the average (of correlation functions) over a finite set of random vectors converges (when we make the size of the set larger and larger) to the value obtained for infinite set. Is it possible to speed this convergence by a suitable choice of sampling algorithm? (may be I need to use some analogue of Metropolis algorithm for discrete distribution).
You could try a low-discrepancy sequence, e.g. with Sobol.jl, instead of random numbers. You could also try one of the importance-sampled Monte–Carlo methods in Cuba.jl.
But essentially you are doing very high-dimensional integration (I’m assuming `N` is large?), and there are no universal solutions for this — high-dimensional integration is intrinsically hard to do accurately via brute-force methods. There are a lot of specialized techniques that people have used on specific problems, but it requires some problem-specific analysis. For example, you say your problem involves modeling disorder — if the disorder is small, and you can use a perturbation-theory technique to linearize the problem as a function of the disorder, then there are a lot of analytical techniques you can employ to simplify matters (often by many orders of magnitude).