mathjax

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$.

No comments:

Post a Comment