How do I simulate flip of biased coin?

Question

In unbiased coin flip H or T occurs 50% of times.

But I want to simulate coin which gives H with probability 'p' and T with probability '(1-p)'.

something like this:

def flip(p):
   '''this function return H with probability p'''
   # do something
   return result

>> [flip(0.8) for i in xrange(10)]
[H,H,T,H,H,H,T,H,H,H]

Solution 1

random.random() returns a uniformly distributed pseudo-random floating point number in the range [0, 1). This number is less than a given number p in the range [0,1) with probability p. Thus:

def flip(p):
    return 'H' if random.random() < p else 'T'

Some experiments:

>>> N = 100
>>> flips = [flip(0.2) for i in xrange(N)]
>>> float(flips.count('H'))/N
0.17999999999999999  # Approximately 20% of the coins are heads

>>> N = 10000
>>> flips = [flip(0.2) for i in xrange(N)]
>>> float(flips.count('H'))/N
0.20549999999999999  # Better approximation 

Solution 2

Do you want the "bias" to be based in symmetric distribuition? Or maybe exponential distribution? Gaussian anyone?

Well, here are all the methods, extracted from random documentation itself.

First, an example of triangular distribution:

print random.triangular(0, 1, 0.7)

random.triangular(low, high, mode):

Return a random floating point number N such that low <= N < high and with the specified mode between those bounds. The low and high bounds default to zero and one. The mode argument defaults to the midpoint between the bounds, giving a symmetric distribution.

random.betavariate(alpha, beta):

Beta distribution. Conditions on the parameters are alpha > 0 and beta > 0. Returned values range between 0 and 1.

random.expovariate(lambd):

Exponential distribution. lambd is 1.0 divided by the desired mean. It should be nonzero. (The parameter would be called “lambda”, but that is a reserved word in Python.) Returned values range from 0 to positive infinity if lambd is positive, and from negative infinity to 0 if lambd is negative.

random.gammavariate(alpha, beta):

Gamma distribution. (Not the gamma function!) Conditions on the parameters are alpha > 0 and beta > 0.

random.gauss(mu, sigma):

Gaussian distribution. mu is the mean, and sigma is the standard deviation. This is slightly faster than the normalvariate() function defined below.

random.lognormvariate(mu, sigma):

Log normal distribution. If you take the natural logarithm of this distribution, you’ll get a normal distribution with mean mu and standard deviation sigma. mu can have any value, and sigma must be greater than zero.

random.normalvariate(mu, sigma):

Normal distribution. mu is the mean, and sigma is the standard deviation.

random.vonmisesvariate(mu, kappa):

mu is the mean angle, expressed in radians between 0 and 2*pi, and kappa is the concentration parameter, which must be greater than or equal to zero. If kappa is equal to zero, this distribution reduces to a uniform random angle over the range 0 to 2*pi.

random.paretovariate(alpha):

Pareto distribution. alpha is the shape parameter.

random.weibullvariate(alpha, beta)

Weibull distribution. alpha is the scale parameter and beta is the shape parameter.

Solution 3

import random
def flip(p):
    return (random.random() < p)

That returns a boolean which you can then use to choose H or T (or choose between any two values) you want. You could also include the choice in the method:

def flip(p):
    if random.random() < p:
        return 'H'
    else:
        return 'T'

but it'd be less generally useful that way.

Solution 4

How about:

import numpy as np
n, p = 1, .33  # n = coins flipped, p = prob of success
s = np.random.binomial(n, p, 100)

Solution 5

One can sample from the X ~ Bernoulli(p) distribution nsamples times using sympy too:

from sympy.stats import Bernoulli, sample_iter
list(sample_iter(Bernoulli('X', 0.8), numsamples=10)) # p = 0.8 and nsamples=10
# [1, 1, 0, 1, 1, 0, 1, 1, 1, 1]

Return 'H' or 'T' instead using

def flip(p, n):
    return list(map(lambda x: 'H' if x==1 else 'T', sample_iter(Bernoulli('X', p), numsamples=n)))

print(flip(0.8, 10)) # p = 0.8 and nsamples=10
# ['H', 'H', 'T', 'H', 'H', 'T', 'H', 'H', 'H', 'H']

View more solutions