mathjax

Showing posts with label Bernoulli. Show all posts
Showing posts with label Bernoulli. Show all posts

Wednesday, January 14, 2015

Bernoulli Distribution

The Bernoulli Distribution models the likelihood of two outcomes: either True or False, Heads or Tails, 1 or 0, Up or Down, Success or Failure, etc. The distribution takes one parameter, $p$, with $0 \leq p \leq 1$, which represents the likelihood of success. Since there are only two outcomes in the sample, $k \in \{0,1\}$, the likelihood of failure is represented by $q=1-p$. This way, the probability of all possible outcomes sums to one: $p+q=1$, and hence it is a legitimate probability distribution.
Formally, if $X$ is a random variable modeled by a Bernoulli distribution with parameter $p$, we say that $X \sim Bern(p)$, $P(X=1)=p$ and that $P(X=0)=1-p=q$.
A Bernoulli experiment or trial involves sampling a Bernoulli distribution once. For example, the flipping of a coin is a Bernoulli trial. If the coin is fair, then $p=0.5$, and by the definition, $q=1-p=0.5$; that is, the chance of heads is the same as the chance of tails.
The $pmf$ is
$p(k;p)=
\left \{
 \begin{array} {c c}
  p & \text{if} \: k=1, \\
  1-p=q & \text{if} \: k=0.
 \end{array}
\right \}
$
The way to read $p(k;p)$ is the probability of $k$, conditioned on $p$. If you look at the Bernoulli distribution as a special case of the Binomial distribution, then $p(k;p)=\binom{1}{1}p^{k}q^{1-k}=p^{k}q^{1-k}$ for $k \in \{0,1\}$.
The expected value of $X$ is $E(X)=p$ and the variance is $Var(X)=p(1-p)=pq$.
import numpy, scipy
import pandas as pd
import matplotlib.pyplot as plt

def bernoulli_trial(p):
    if numpy.random.random() < p:
        return 1
    return 0

p = 0.1
N = 20
trials = [bernoulli_trial(p) for i in range(N)]
print(trials)
This results in the following list:
[0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
Changing $p$ to $0.5$ gives something like:
[0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1]
Plot like so:
plt.scatter(range(len(trials)), trials, linewidth=0, alpha=0.5, s=50)
plt.title('Bernoulli Trials')
plt.xlabel('trial number')
plt.ylabel('trial outcome')
plt.show()
This data is called a random sample. We can use maximum likelihood estimation to estimate the value of $p$ if we don't know it.
The next piece of code does 500 Bernoulli trials for 9 different values of $p$: $p \in \{0.1, 0.2, ..., 0.9 \}$. We will pretend we don't know $p$ for each of these 9 samples.  
plt.clf()
N = 500
data = []
for p in [float(i)/10.0 for i in range(1,10)]:
    trials = [bernoulli_trial(p) for i in range(N)]
    data.append(trials)
data = numpy.array(data).T
df = pd.DataFrame(data, columns = [float(i)/10.0 for i in range(1,10)])
means = df.mean()
means.plot(kind='bar')
plt.xlabel('actual $p$')
plt.ylabel('observed $p$')


As you can see, the observed $p$ closely approximates the actual $p$.