Looking for a function that yields a value with a certain probability?

Questions
1

I need a function that returns True or 1 or whatever for a given probability. For probabilities of the form p=1/n where n is an integer this is easy but for p=m/n this seems to be more difficult. For p=3/10 it could be achieved with

r=RandomInteger[{1, 10}];
If[r==1||r==2||r==3,c=1,c=0]

So c=1 will be returned with p=3/10 but I am wondering if there is not a standard function for that purpose. Any ideas?

2

In[7]:= f[prob_, n_Integer, seed_Integer] :=
Block[{}, SeedRandom[seed];
Boole[# <= prob] & /@ RandomReal[{0, 1}, n]]

In[9]:= f[.3, 20, 0]

Out[9]= {0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0}

In[10]:= f[.3, 20, 1]

Out[10]= {0, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1}
3

Here are two ways to get a list of 0’s and 1’s with a specified probability:

p = 0.375;  (* Probability of observing a 1 *)
n = 10;  (* Number of random samples *)

x = RandomChoice[{p, 1 - p} -> {1, 0}, n]
(* {0,0,0,0,1,0,1,1,0,0} *)

x = RandomVariate[BernoulliDistribution[p], n]
(* {0,0,1,0,1,1,1,0,0,0} *)
4

Thanks a lot for your help. I always learn something new from your proposals. I tested your methods for speed and finally found that the most rapid code for returning c=1 with probability p=m/n is this one:

   If[1 <= RandomInteger[{1, n}] <= m, c = 1, c = 0];
5

RandomChoice is indeed the way to go!

6

Would you elaborate on how you tested the speed? I’m assuming it must be in how you generated multiple samples from

If[1 <= RandomInteger[{1, n}] <= m, c = 1, c = 0];

I get the following:

nSimulations = 10000000;
n = 10;
m = 7;
p = m/n;(*Probability of observing a 1 *)

Timing[x = Table[If[1 <= RandomInteger[{1, n}] <= m, c = 1, c = 0], {i, nSimulations}];]
(* {10.140625`,Null} *)

Timing[x = RandomChoice[{p, 1 - p} -> {1, 0}, nSimulations];]
(* {0.171875`,Null} *)

Timing[x = RandomVariate[BernoulliDistribution[p], nSimulations];]
(* {0.171875`,Null} *)
7

I tested the function with this code as I use it in a Do loop similar to this one:

z = 10^6; c = 0;
Do[
  If[1 <= RandomInteger[{1, 10}] <= 3, c++],
  {z}] // Timing

Output is {0.920406, Null}

z = 10^6; c = 0; p = 3/10;
Do[
  If[RandomChoice[{p, 1 - p} -> {1, 0}] == 1, c++],
  {z}] // Timing

Output is {4.43043, Null}

Of course, in both cases c/z tends to 0.3.

8

If you just need to total number of 1’s and don’t need to keep the individual observations, you might consider the following which is a lot faster:

Total[RandomChoice[{p, 1 - p} -> {1, 0}, nSimulations]] // Timing
(* {0.046875, 699861} *)

And in that case even faster is the following:

RandomVariate[BinomialDistribution[z, p], 1] // Timing
(* {0., {700138}} *)