MATLAB: How to similuate a coin flip with probablility p

probability

How do I simulate getting a result, either 0 or 1, with probability p. So if p=0.5 I should get an output of 0 half of the time, and 1 half of the time.

Best Answer

100 tosses with p=0.5.

x = round(rand(100,1));

If you want a probability other than p=0.5, then realize that rand() is uniform random number generator between [0,1], so you can assign the output of rand() accordingly. For example, for p=0.25:

y = zeros(100,1);
x = rand(100,1);
y(x<0.25) = 1;

Or if you just want to simulate the number of 0's or 1's in a certain number of trials. Say 100 for example. Here is a simulation of ten such experiments. Requires Statistics Toolbox.

R = binornd(100,0.5,10,1);

Related Solutions

MATLAB: Simulating an unbias coin with a bias coin flip

I think that the simplest approach is to write a little function to simulate a single toss, and then build up your simulation from that. Here is that function:

function side = simulateOneToss
% Make two tosses (outcomes 0 and 1 could stand for heads and tails)

twoTosses = round(rand(1,2));
 % If outcome is not HT or TH (both of which sum to 1), try again.

while sum(twoTosses) ~= 1
    twoTosses = round(rand(1,2));
end
 % Take first of the two tosses as the answer.

side = twoTosses(1);

Based on the comments below, here is a version that allows a probability other than 0.5:

function side = simulateOneToss(p)
% Make two tosses (outcomes 0 and 1 could stand for heads and tails)
twoTosses = rand(1,2)<p;
% If outcome is not HT or TH (both of which sum to 1), try again.
while sum(twoTosses) ~= 1
  twoTosses = rand(1,2)<p;
end
% Take first of the two tosses as the answer.
side = twoTosses(1);

Note that I don't round before comparing with the probability (otherwise you'll always get p=0.5 in effect). Now this can be used in a simulator:

function outcomes = simulateTosses(p,nTosses,nTrials)
tosses = zeros(nTosses,nTrials);
for i=1:numel(tosses)
  tosses(i) = simulateOneToss(p);
end
outcomes = sum(tosses>p);

Finally, you can run it and plot the results as follows:

nTosses = 100;
nTrials = 10;
p = 0.6;
outcomes = simulateTosses(p,nTosses,nTrials);
bar(1:nTrials,outcomes)

MATLAB: How to interpret the outputs of the function “runstest.m”

This is a problem of appreciating the terminology, and the double negative way things get worded. Do two negatives make a positive? Not always. But then, not, not, always either. ;-)

Lets see what happens when we throw a very non-random sequence at runstest.

runstest(cumsum(rand(1,100)))
ans =
     1
runstest(mod(1:100,2))
ans =
     1

The null hypothesis is that the sequence posed is truly random. Can we KNOW a sequence was a random one? Well, you can't KNOW that as fact. You can only test how random it seems to be, or, how non-random. So both cases above are clearly seen to be unlikely events, IF the sequence were truly a random one. Essentially, we reject the null hypothesis, so the a priori assumption that the sequences really are random.

But, when I do this:

runstest(rand(1,1000))
ans =
     0

runstest is not willing to say the sequence seems to be non-random. So, we cannot reject the null hypothesis. Does that mean the sequence is random? NO!!!!!! Only that we cannot comfortably decide the sequence is probably non-random.

So, now lets look at p.

[h,p] = runstest(cumsum(rand(1,100)))
h =
     1
p =
      4.04395443273876e-29

Yep. It is non-random. I'll bet some decent money on that. p is TINY, it reflects a measure of how likely we think this event might be, in context of a truly random sequence.

Now, lets try a sequence of sequences.

x = rand(1,100);
[h,p] = runstest(x)
h =
     0
p =
         0.838684318956303
[h,p] = runstest([x,x])
h =
     0
p =
         0.720356034109542
[h,p] = runstest([x,x,x])
h =
     0
p =
          0.64038079177353
[h,p] = runstest([x,x,x,x])
h =
     0
p =
         0.578177151383335

See that the way I've constructed those sequences, the latter ones are seen to be somewhat less likely to be random in context of the runs seen, although runstest is not designed to detect the totally non-random way I extended the sequences. It cannot detect long lags as I extended the sequence, which in fact, make the longer versions non-random.

If I try a sequence of length 400 directly from rand, p is now larger.

[h,p] = runstest(rand(1,400))
h =
     0
p =
          0.72633274852919

So h==1 says runstest has detected an event that did not seem consistent with randomness. h==0 says it did not detect significant non-randomness, but it cannot tell you to conclude a sequence is truly random, or not. p tells you the strength of what it sees.

Best Answer

Related Solutions

MATLAB: Simulating an unbias coin with a bias coin flip

MATLAB: How to interpret the outputs of the function “runstest.m”

Related Question