Basic Concepts in Probability and Information Theory

Contents

1. References
2. Probability Theory
   1. Random Variables
   2. Independence
   3. Continuous case: Probability densities
   4. Expectation and covariance
3. What can be done with a distribution
4. Likelihood and Maximum Likelihood Estimation
5. Useful Distributions
   1. Bernouilli binary distribution
   2. Binomial Distribution
   3. Choosing a Prior for the Binomial Distribution
   4. Deriving the MAP Estimate of the Binomial Distribution with a Beta Prior
   5. Bayesian Estimate of the Binomial
   6. Multinomial Distribution
   7. Posterior Predictive Distribution
   8. Bayesian Method: Model, Observe, Update, Predict
6. Decision Theory
   1. MAP estimate minimizes 0-1 loss
   2. Reject Option
7. Central Limit Theorem
8. Information Theory
   1. Entropy
   2. Entropy and Huffman Coding
   3. Cross Entropy and KL Divergence

References

1. Chapter 1 of Pattern Recognition and Machine Learning by Chris Bishop
2. Machine Learning: a Probabilistic Perspective, by Kevin Murphy 2012
3. Chapter 2 of Machine Learning by Tom Mitchell, 2016, on "Estimating Probabilities" provides a detailed view of the computation of MLE and MAP for the Beta-Binomial model discussed below (available online).

1. Probability Theory
2. Random Variables
3. Information Theory
4. Statistical Inference
5. Markov Models
6. Parts of Speech Tagging: Statistical Models

In particular, read:

1. Probabilistic and Bayesian Analytics
2. Cross-validation for detecting and preventing overfitting

Probability Theory

Consider the set Ω of outcomes. An event is defined as a subset of outcomes from Ω -- there are 2ⁿ events over Ω. Define a function p that meets the following conditions:

p: 2Ω → [0..1] (that is, for all event e, 0 ≤ p(e) ≤ 1)
p(A ∪ B) = p(A) + p(B) - p(A ∩ B)
p(∅) = 0
p(Ω) = 1

Random Variables

The computations predicted by probabilities become interesting and sometimes surprising when we consider complex experiments, involving multiple random variables. The dependency relations among the variables are reflected in the computations using simple rules -- called "sum and product rules of probability".

Consider the classical example of boxes and balls with Bishop's figure:

In this setting, we are considering a "structured experiment" which proceeds in 2 stages: first pick a box (either red or blue), then pick a ball within the selected box (either orange or green). Assume we repeat this experiment many times, and assume we get the following aggregated observations:

1. p(red box) = 40%
2. p(blue box) = 60%
3. All balls in the box are equi-probable (can be picked equally likely)

1. B ∈ {blue, red} (the box can be either blue or red)
2. F ∈ {orange, green} (the ball can be either orange or green)

In this view, we can define several interesting quantities (assuming we ran the experiment N times, and N → ∞):

• p(X = x_i, Y = y_j) = n_ij/N is the joint probability
• p(X = x_i) = Σ_j p(X = x_i, Y = y_j) = c_i/N is the marginal probability
• p(X = x_i | Y = y_j) = n_ij / y_j is the conditional probability (read as probability of X=x_i given that Y=y_j)

p(X = xi, Y = yj) = p(Y = yj | X = xi)p(X = xi)
p(X = xi, Y = yj) = p(X = xi | Y = yj)p(Y = yj)

p(Y | X) = p(X | Y) p(Y) / p(X)
p(X) = ΣYp(X|Y)p(Y)

p(B = r) = 4/10
p(B = b) = 6/10

p(F = g | B = r) = 1/4
p(F = o | B = r) = 3/4

p(F = g | B = b) = 3/4
p(F = o | B = b) = 1/4

p(F = o) = p(F = o | B = r)p(B = r) + p(F = o | B = b)p(B = b) = 9/20 = 0.45
p(F = g) = 11/20 = 0.55

Another non-trivial computation is obtained by applying Bayes rule. Bayes rule in general is used to "reverse dependency". Consider the following observation: we run the experiment and we obtain an orange ball. What is the likelihood that the ball was picked from the red box?

p(B = r | F = o) = p(F = o | B = r)p(B = r)/p(F = o) = 2/3

You can solve the example given in the article quoted above: An Intuitive Explanation of Bayes' Theorem using this formula: 1% of women at age forty who participate in routine screening have breast cancer. 80% of women with breast cancer will get positive mammographies. 9.6% of women without breast cancer will also get positive mammographies. A woman in this age group had a positive mammography in a routine screening. What is the probability that she actually has breast cancer? The solution is 7.76% -- which is quite not intuitive even on a very simple binary case. Make sure to work out the formula to convince yourself of how the result is obtained.

Independence

p(X,Y) = p(X)p(Y)

In general:
p(X, Y) = p(X | Y) p(Y)

When X and Y are independent, then:
p(X, Y) = p(X | Y) p(Y) = p(X) p(Y)

Hence:
p(X | Y) = p(X)

We can extend the independence relation to conditional independence:

X and Y are conditionally independent in Z iff:
p(X,Y|Z) = p(X|Z)p(Y|Z)

Continuous case: Probability densities

p(x ∈ [a..b]) = ∫ab p(x)dx

∀ x, p(x) ≥ 0
∫-∞+∞ p(x)dx = 1

p(x) = ∫ p(x,y)dy
p(x,y) = p(y|x)p(x)

Expectation and covariance

Ex(f) = Σx p(x)f(x)  for the discrete case
Ex(f) = ∫x p(x)f(x)dx  for the continuous case

If we sample x many times and obtain values (x₁...x_N), then the following approximation holds:

E(f) ≅ 1/N Σn=1Nf(xn)

Ex[f(x,y)] = Σxf(x,y)  (This is a function of y)
Ex[f | y] = Σx p(x | y) f(x)

var(f) = E[(f(x) - E(f(x)))2] = E[f(x)2] - E[f(x)]2

cov[x,y] = Ex,y[{x - E(x)}{y - E(y)}] = Ex,y[xy] - E[x]E[y]

Expectation and variance are useful summaries of a complex probability distribution -- they tell us in 2 numbers something useful about a full distribution -- what is the mean value and how much variability there is around the mean. Often, we want to model the following situation: we know that a quantity has been measured with some level of noise; in such a case, we are looking for a way to describe what we know about the quantity. One way is to state that the actual quantity has a certain expected value (derived from the measurement) with some level of variance (corresponding to the accuracy of the measurement method). In this case, the 2 quantities that characterize the possible distribution of the quantity are its expectation and its variance. If we know more about the measurement process, we can suggest a finer distribution model. Otherwise, we use a standard distribution model (typically, we use the Gaussian distribution in such cases as described below).

In the case of two variables, covariance is a measure of how much two random variables change together. For example, if in general we observe that when x grows, then y tends to grow as well, and when x decreases, then y decreases as well, that is, x and y show in general a similar pattern, then the covariance is positive. If the behavior is inverse, that is, when x grows, y decreases and vice versa, then the covariance is negative. If there is no correlation between x and y (they are independent), then the covariance is 0. The interpretation of the magnitude of the covariance is not easily interpreted.

What can be done with a distribution

The following are the key operations:

• Sample: draw a random value from the distribution, such that, if we keep sampling many times, the sampled objects are distributed as predicted by the distribution.
• Estimate: given a set of samples, guess which distribution is likely to have generated these samples. Estimation is the reverse operation of sampling. We can think of it as a "factory method".
• Marginalize: given a joint distribution p(x,y), compute the marginal distribution p(x) or p(y)
• Condition: given a joint distribution p(x, y), compute the conditional distribution p(x|y) or p(y|x)
• Expectation: compute the expectation of the distribution Ep(x) - think of it as the weighted average of the distribution
• Variance: compute the variance of the distribution - think of it as the measure of variation around the average
• Average: given a function f, compute the weighted average of the function over the distribution, that is Ep(f(x))

When X is drawn from a distribution D, we note:

X ~ D

In Python, the scipy module provides an implementation and an API to manipulate distributions - both discrete and continuous. See the tutorial on the scipy.stats module. See the Scipy API for rv_discrete and rv_continuous - which are the two abstract interfaces that define Distributions in Scipy.

This notebook on Scipy Discrete Distributions (ipynb) illustrates how to use such an interface to construct the distribution object, sample from it, compute the expected value for various functions and the variance and standard deviation.

Likelihood and Maximum Likelihood Estimation

p(w|D) = p(D|w)p(w)/p(D)

In the Bayes equation, we use the following terminology:

• Posterior: p(w|D) is the updated distribution over w given the observation D (that is, it reflects what we know about w after we have observed D).
• Prior: p(w) is the distribution over w independently of observations (that is, it reflects what we know about the parameters independently of any observation).
• Likelihood: p(D|w) is considered as a function over w, it reflects how likely it is that we sampled the dataset D from the distribution given a value of w. Note that the likelihood as a function of w is NOT a distribution over w (that is, there is no normalization constraint that the sum of the likelihood for all possible values of w be 1).

posterior ∝ likelihood x prior

p(D) = ∫ p(D|w)p(w)dw

The frequentist approach to estimation consists in finding the best possible value for w given an observation D. In this case, the normalization constant p(D) plays no role, and the output is an optimal value for w, w*. That is, w* is a pointwise estimate of the distribution parameter. The Bayesian approach instead tries to estimate the distribution over w, p(w|D). This requires the computation of the normalization constant p(D).

In general, Bayesian approaches to estimation are most useful if the task we perform (regression or classification) is one step within a multi-step process. The Bayesian approach is preferable in this case because the output of the classification for example, does not need to choose the most likely value for w, but instead can let further steps choose among multiple options that may be likely and can be determined based on features which are not considered at the time of the classification. For example, consider the task of syntactic parsing. Parsing often uses as features the parts of speech of each word. Therefore, the classification task of parts of speech tagging is one step in the 2-step process of parsing. The frequentist approach will compute the most likely POS tag for each word before parsing starts. The Bayesian approach will instead compute a distribution over possible POS tags for each word, which can be later used by the parser.

Estimation can proceed in 3 different ways, of increasing complexity:

1. Maximum likelihood: find the parameter w* which yields the maximum likelihood. w* = argmax p(D|w). This ignores the prior term p(w) and the normalization constant p(D) in Bayes's equation. Maximum likelihood estimators have a tendency to over-fit the data.
2. Maximum a posteriori estimate (MAP): this method tries to maximize the likelihood together with a prior: find the parameter w* which yields the maximum value for the composed term likelihood x prior, assuming a given prior distribution over w. w* = argmax p(D|w)p(w). This has the effect of avoiding overfitting and is similar to regularization. This ignores the normalization constant p(D) in the Bayes equation.
3. Bayesian estimation: find the posterior distribution p(w|D) by computing the full right-hand side of Bayes equation. p(w|D) = p(D|w)p(w)/∫ p(D|w)p(w)dw. This is the most informative way to estimate, but is computationally expensive.

Useful Distributions

When solving a supervised machine learning problem, we attempt to characterize the distribution of a random variable (the value to predict) given a set of examples (observations). One way to characterize this operation is as an operation of estimation of the parameters of a distribution from a known family given the observations. We will, therefore, often end up performing continuous parameter estimation (and use integrals) to learn a function that maps input parameters to predicted values on the basis of discrete examples.

Bernouilli binary distribution

p(x=1|μ)=μ  where 0 ≤ μ ≤ 1
p(x=0|μ)=1-μ

Bern(x|μ)=μx(1-μ)(1-x)

E(p) = μ
var(p) = μ(1-μ)

Computing the maximum likelihood estimate of a Bernouilli distribution

Given D = (x1...xN):

p(D|μ) = ∏n=1Np(xn|μ) = ∏ μxn(1-μ)(1-xn)

We derive this expression from the fact that the observations are 
independent of each other (hence the product).

In order to derive the maximum likelihood estimate μML, we must solve:
μML = argmaxμ p(D|μ)

To obtain the maximum of the likelihood function, we must compute its derivation.  
It is unpleasant to derive long products, so instead, we will look for the maximum 
of the log-likelihood function. The log likelihood is:

ln p(D|μ) = Σ ln p(xn|μ) = Σ {xnlnμ + (1-xn)ln(1-μ)}

We take advantage of two properties of the log function here:

Log is a monotonic function - and therefore, argmax(log(f)) = argmax(f)
Log(axb) = Log(a) + Log(b) (and hence Log(an) = nLog(a))

Remember also that d ln(x)/dx = 1/x.


By computing the derivation of the log likelihood as a function of μ,
we derive that the likelihood reaches its maximum on:

d LogLikelihood(μ) / dμ = Σ {xn/μ - (1-xn)/(1-μ)}

μML = 1/N Σ xn

With a dataset D=(x₁...x_N), we assume that for all i in (1..N), x_i ~ Bernouilli(μ) and we say that the N variables x_i are independent and identically distributed (IID). As as typical for MLE (maximum likelihood estimates), this MLE estimate severely overfits the distribution when the dataset is small.

Binomial Distribution

Bin(m|N,μ) ∝ μm(1-μ)N-m

Bin(m|N,μ) = (N m) μm(1-μ)N-m

E[m] = Nμ
var[m] = Nμ(1-μ)

Choosing a Prior for the Binomial Distribution

To this end, note that the likelihood function of the Bernouilli distribution has a form which is proportional to μ^x(1-μ)^1-x. If we choose the prior to also be a power of μ and (1-μ), then the posterior distribution will also have the same form as the prior. This is convenient because it means the posterior distribution will belong to the same family of functions as the prior, and the Bayesian update procedure will consist of only updating the parameters of the prior distribution. This convenient property is called conjugacy: we identify the conjugate prior of the binomial distribution to be:

Beta(μ|a,b) = Γ(a+b) / Γ(a)Γ(b) μa-1(1-μ)b-1

where Γ(x) = ∫ ux-1e-udu

E[μ] = a / (a+b)
var[μ] = ab / (a+b)2(a+b+1)

Deriving the MAP Estimate of the Binomial Distribution with a Beta Prior

μMAP = argmax μ [p(μ|D)] = argmax μ [p(D|μ)p(μ)]
μMAP = argmax μ [log p(D|μ) + log p(μ)]

Recall that:
μMLE = argmax μ [log p(D|μ)]

Let us now compute μ_MAP with a Beta prior:

p(μ) = Beta(μ | a, b) = Γ(a+b) / Γ(a)Γ(b) μa-1(1-μ)b-1
Bin(m|N,μ) = (N m) μm(1-μ)N-m

μ ~ Beta(a, b)
m ~ Bin(μ, N)

μMAP = argmax μ [log p(D|μ) + log p(μ)]

We compute the derivation of the terms:
l(μ) = p(D|μ) * p(μ) = (N m) μm(1-μ)N-m * Γ(a+b) / Γ(a)Γ(b) μa-1(1-μ)b-1
l(μ) = (N m) Γ(a+b) / Γ(a)Γ(b) μm+a-1(1-μ)N-m+b-1

log l(μ) = log [(N m) Γ(a+b) / Γ(a)Γ(b)] + (m+a-1)log(μ) + (N-m+b-1)log(1-μ)

d log l(μ) / dμ = (m+a-1)/μ - (N-m+b-1)/(1-μ)

μMAP = (m + a - 1) / (N + a + b - 2)

Example: Assume we toss a coin 3 times and obtain 3 times Tails. In this case, D = (0,0,0).

μ_MLE = 0 (a very extreme conclusion -- the coin will always drop on Tails).

If we take a prior Beta(2,3) - see the bottom left graph above, which indicates a preference with a peak preference for μ = 0.4 -- then μ_MAP = 1/6 (a less extreme conclusion).

If we take a prior Beta(8,4) - see the bottom right graph above, which indicates a preference with a peak preference for μ = 0.6 -- then μ_MAP = 7/13, our peak preference is reduced, but remains quite high, because our prior was very strong, and 3 observations were not sufficient to force our conclusion more.

Bayesian Estimate of the Binomial

If we choose a Beta conjugate prior for μ, we derive from Bayes theorem:

Assume:
Prior μ ~ Beta(a, b)
      X ~ Bin(μ, N)

p(μ|N0,N1,a,b) ∝ μN1+a-1(1-μ)N0+b-1

where N1 is the number of 1s and N0 = N-N1 is the number of 0s observed in the dataset D,
and a and b are hyper-parameters that determine our prior knowledge about p(μ).

The posterior of the binomial given a dataset with (m draws of 1, z=N-m draws of 0)
and a prior Beta(a,b) is: Beta(a+m, b+z)

Another way of saying this is to imagine a sequential Bayesian update process: assume that instead of updating our posterior belief for a whole dataset of N observations at once, we update our belief sequentially, for each observation point, one by one. In this way, our prior can start from a uniform distribution (which happens to be Beta(1,1)), and for each datapoint, we update the distribution by adding one to either the a hyper-parameter (if the datapoint is 1) or to the b hyper-parameter (if the datapoint is 0).

This process of sequential update is illustrated graphically below: The graph on the left shows the prior distribution over the μ parameter: it indicates a "rough" preference for a fair Bernouilli experiment (E(μ)=0.5) with high variance (low certainty) encoded as a Beta(2,2) prior distribution. The likelihood function is shown in the center diagram, corresponds to a single observation of N=N₁=1 (in which case, the MLE estimate is μ = 1). The posterior updates the prior distribution from Beta(2,2) to Beta(3,2), which shifts smoothly from 0.5 towards 1, and has lower variance than the original Beta(2,2) -- but still does not "jump" to the extreme overfitted conclusion that μ=1.

Remember that the mean of Beta(a,b) is a/(a+b) and the variance is ab/(a+b)²(a+b+1). Given these we can encode our prior beliefs in the following manner: assume we want a prior for our distribution that has a mean of 0.7 and a standard deviation of 0.2. We resolve for a,b in the equations:

mean = 0.7 = ab/(a+b)
var = 0.22 = ab/(a+b)2(a+b+1)

The following Notebook on Bernouilli Bayesian Estimate (code) illustrates how to perform such bayesian estimation using Scipy code.

Multinomial Distribution

v = (0 0 ... 0 1 0 ... 0): vi = 1 and vj=0 for all j≠i for i,j ∈ (1..K)
encodes the outcome i out of the 1 to K possible outcomes.

p(x | μ) = ∏μkxk

where μ = (μ1, ..., μK) is the parameter of the distribution and

Σμk = 1

E(x|μ) = Σp(x|μ)x = μ

Mult(m1, m2, .., mK | N, μ) = (N choose  m1 ... mK)∏μkmk
      = N! / m1! ... mK!∏μkmk

Consider the way we can compute the likelihood of a dataset D = (x₁ ... x_n) where each x_i = (x_i1 ... x_iK) is the outcome of a trial. Then the likelihood of the dataset D is given as follows:

likelihoodD(μ) = p(D|μ) = ∏n=1..N ∏k=1..K μkxnk = ∏k=1..K μΣxnk

likelihoodD(μ) = p(D|μ) = ∏k=1..K μkmk

log-likelihood(μ) = Σk=1..K mk ln(μk)

For Bayesian inference, we introduce a conjugate prior on the multinomial which is a direct generalization of the Beta distribution introduced for the binomial distribution. This conjugate distribution is called the Dirichlet distribution and is defined as:

Dirichlet(μ | α) = Γ(α0) / Γ(α1) ... Γ(αK) ∏μkαk-1

where α0 = Σαk
and   μ = (μ1,..., μK)

p(μ | D, α) ∝ p(D | μ) p(μ | α) ∝ ∏μkαk+mk-1

and hence, exploiting the conjugacy of the Dirichlet prior:

p(μ | D, α) = Dirichlet(μ | α+m)

Posterior Predictive Distribution

In a Bayesian context, the posterior predictive distribution makes use of the entire posterior distribution of the parameters given the observed data to yield a probability distribution over an interval rather than simply a point estimate. It is computed by marginalising over the parameters, using the posterior distribution:

p(d | D, α) = ∫w p(d | w) p(w | D, α) dw

p(d | D, α) = Ew| D, α[ p(d|w) ]

Bayesian Method: Model, Observe, Update, Predict

The method used to estimate the joint distribution is supervised:

1. We observe a dataset of labelled (features, label) pairs.
2. We hypothesize that the distribution follows a probabilistic model with a set of parameters - that is, we hypothesize that the observed pairs have been sampled from a specific distribution model with specific parameters.
3. We encode knowledge about the process that has produced the observed data in the form of a prior distribution over the probabilistic model parameters.
4. We estimate the parameters by using one of the estimation methods discussed above (MLE, MAP, Posterior distribution over the parameters).
5. We then derive a new posterior predictive distribution that takes into account the model, the prior and the observations.
6. On the basis of the posterior predictive distribution, we can assess p(features, label) for any new unseen features vector and any candidate label.

Example: Language Model using Bag of Words (Murphy 2012, p81)

Remember the spell checker example discussed in the first lecture.
We discussed there the noisy channel model over words and modelled spell checking as a probabilistic model:

argmaxc P(w|c) P(c)

That is: given an observed word w and a set of candidate words {c...} which are the possible words intended by the writer,
the spelling correction of w is the candidate word c which maximizes the product p(w|c)*p(c).

In this term - we called:
p(c): the language model - how likely it is that the word c occurs in a sentence
p(w|c): the error model - how likely it is that w was written when the writer intended c.

We will design here a Bayesian model to estimate the distribution p(c) - the language model, using the Dirichlet-multinomial model.

In general, a language model predicts which words occur in a sequence of words (such as a document or a sentence).
We model each word is a choice of a word within a vocabulary i in {1...K}.
In the simplest model - we consider that the i-th word in the sequence is sampled independently from all the other words:

p(Xi | X1...Xi-1) = p(Xi)

which we write as:

Xi ~ Cat(w)

Where Cat(w) is a categorical distribution over K possibilities - which is the special case of 
a multinomial distribution with a single draw.

Cat(X | w) = Multinomial(X | 1, w) = Πj=1..KwjI(xj=1) 
where I(x = 1) is the indicator function which is 1 when x=1, 0 otherwise.

We observe a dataset - that is a sequence of words: 

Mary had a little lamb, little lamb, little lamb,
Mary had a little lamb, its fleece as white as snow


We first represent the vocabulary as K-dimensional vectors with one-hot encoding:

Word	Index
=================
mary	1
lamb	2
little	3
big	4
fleece	5
white	6
black	7
snow	8
rain	9
unk	10

The token unk is used to encode any "unknown" word that may appear at prediction or further training time.

We then encode the observed dataset: we remove punctuation, and very common words ("a", "as", "the" etc).


1 10 3 2 3 2 3 2
1 10 3 2 10 5 10 6 8


We now ignore the word order and summarize this dataset by simply counting the number of occurrences of 
each word in the vocabulary:

Word	Index	Count
=========================
mary	1		2
lamb	2		4
little	3		4
big	4		0
fleece	5		1
white	6		1
black	7		0
snow	8		1
rain	9		0
unk	10		4

This is a histogram of word counts in the observed dataset.  This model of text is called bag of words
because it ignores the order in which words occur in the text.  The dataset is overall encoded as a vector
of counts D = (N1, ..., N10).

Our objective is now to estimate the parameters w of the distribution Cat(w) from which this dataset would
have been most likely sampled.  w is also a vector of 10 numbers (w1...w10),

We use a prior Dir(α) for the parameter w - where α is also a vector (α1...α10)
which captures our a prior knowledge about the possible values of the parameter w.

p(X = j | D) = αj + Nj / Σk αk + Nk

If we set α = (1...1), we obtain:

p(X = j | D) = (1 + Nj) / (10 + 17)
             = (3/27, 5/27, 5/27, 1/27, 2/27, 2/27, 1/27, 2/27, 1/27, 5/27)
			 
Observe that the predicted frequency of the words "black" and "big" are non-zero even though they were not observed in the 
dataset.  This is an example of smoothing.

Decision Theory

In the general case, we follow this scheme where we try to "guess" Nature's behavior:

• "Nature" picks a label y in Y - unknown to us.
• Nature generates an observation x in X which we can see.
• We have to make a decision - that is choose an action a from an action space A.
• We then incur some loss L(y,a) which measures how incompatible our action a is with Nature's hidden state y.

Our goal is to design a decision procedure which is a function (δ: X → A) which specifies the optimal action for each possible input. Optimality is measured by the expected loss:

δ(x) = argmina in AE[L(y,a)]

ρ(a | x) = Ep(y|x)[L(y,a)] = ΣyL(y,a)p(y|x)

δ(x) = argmina in Aρ(a|x)

MAP estimate minimizes 0-1 loss

L(y, a) = I(y != a) (1 if y != a, 0 if y = a).

ρ(a|x) = p(a != y | x) = 1 - p(y|x)

y*(x) = argmaxy in Yp(y|x)

Reject Option

The way to model this situation is by defining a loss function:

Assume we have to choose one label a in (1..C).  Define C+1 as the option of "reject" 
(meaning we prefer to pass and do not make a decision).

L(y = j, a = i) = 0 	if i=j and i,j in (1..C)
                  λr	if i = C+1
		  λs	otherwise
				  
Where λr is the cost of rejecting a decision and λs is the cost of a substitution (giving a wrong answer).

Central Limit Theorem

Xi ~ D(μ σ2) are i.i.d from any distribution that has E[D] = μ and var[D] = σ2.

Consider Sn = ΣXi 

Then as n → +∞ p(Sn = x) ~ (2 Π n σ2)-1/2 exp(-(x - Nμ)2/2Nσ2)

In other words:

Zn = (Sn - Nμ)/σN1/2 ~ N(0, 1)

Exercise: Write Python code that draws the density of x ~ Beta(1,5) and of (1/N Σ x) for N=5, 10, 20.

Information Theory

Entropy

H(X) = - Σk p(X = k) log2(p(X=k)

Entropy and Huffman Coding

1. Fixed-length codes: all symbols in the alphabet are encoded with a string of bits of the same length.
2. Variable-length codes: symbols are encoded by codes of different lengths.

Huffman's coding is such a code: it is a variable length code that is also a prefix code: that is, no codeword is a prefix of any other codeword. This property makes it possible to decode the code - the bit encoding of a string is unambiguous.

For example:

A = 0, B=10, C=110, D=111

010011001110 -- 0 | 10 | 0 | 110 | 0 | 111 | 0
             -- A   B    A   C     A   D    A

Consider the following distribution over the symbols A-D:

p(A) = 0.25
p(B) = 0.25
p(C) = 0.20
p(D) = 0.15
p(E) = 0.15

p(B) = 0.25 --
              + 0.45 ---------
p(C) = 0.20 --                |
p(A) = 0.25 ----------|       + 1.0
p(D) = 0.15 --        + 0.55 -|
              + 0.30 -|
p(E) = 0.15 --

B = 00
C = 01
A = 10
D = 110
E = 111

Here is Python code to compute the Huffman code given a distribution over symbols. To follow the code - remember that a heapq (which implements a priority queue) is an array for which heap[k] <= heap[2*k+1] and heap[k] <= heap[2*k+2] for all k, counting elements from zero. The interesting property of a heap is that heap[0] is always its smallest element. heappop(heap) pops out the smallest element out of a queue, and heappush(heap, element) pushes an element into the queue while maintaining the heap invariant.

In the code below, the heap contains items of the form [weight [letter code] ...] where weight is the overall count of all letters that appear in the list, and for each letter, we incrementally maintain the Huffman code of the letter.

# From http://rosettacode.org/wiki/Huffman_codes
from heapq import heappush, heappop, heapify
from collections import defaultdict
 
def encode(symb2freq):
    """Huffman encode the given dict mapping symbols to weights"""
    heap = [[wt, [sym, ""]] for sym, wt in symb2freq.items()]
    heapify(heap)
    while len(heap) > 1:
        lo = heappop(heap)
        hi = heappop(heap)
        for pair in lo[1:]:
            pair[1] = '0' + pair[1]
        for pair in hi[1:]:
            pair[1] = '1' + pair[1]
        heappush(heap, [lo[0] + hi[0]] + lo[1:] + hi[1:])
    return sorted(heappop(heap)[1:], key=lambda p: (len(p[-1]), p))
 
txt = "this is an example for huffman encoding"
symb2freq = defaultdict(int)
for ch in txt:
    symb2freq[ch] += 1
huff = encode(symb2freq)
print("Symbol\tWeight\tHuffman Code")
for p in huff:
    print("%s\t%s\t%s" % (p[0], symb2freq[p[0]], p[1]))

Symbol	Weight	Huffman Code
 	6	101
n	4	010
a	3	1001
e	3	1100
f	3	1101
h	2	0001
i	3	1110
m	2	0010
o	2	0011
s	2	0111
g	1	00000
l	1	00001
p	1	01100
r	1	01101
t	1	10000
u	1	10001
x	1	11110
c	1	111110
d	1	111111

Cross Entropy and KL Divergence

KL(p||q) = Σk p(k) log[p(k) / q(k)]

We can rewrite KL(p||q) as:

KL(p||q) = Σk p(k) logp(k) - Σk p(k) log q(k)
         = E(log p) - E(log q)
		 = - H(p) + H(p, q)

The intuitive meaning of the KL(p||q) divergence can be understood as follows:

• Assume a language generates strings over an alphabet with distribution p (that is, p(a) is the likelihood that symbol a is used in strings of the language).
• We construct an estimate of p that is inaccurate - in the form of the distribution q over the same alphabet.
• We now construct a Huffman encoding over the alphabet on the basis of q.
• The average number of bits per symbol "wasted" because we used q instead of p to encode the symbols over the language is KL(p||q).

Similar distances between distributions exist - most significantly the Hellinger Distance.

"""
Three ways of computing the Hellinger distance between two discrete
probability distributions using NumPy and SciPy.
From https://gist.github.com/larsmans/3116927
"""

import numpy as np
from scipy.linalg import norm
from scipy.spatial.distance import euclidean


_SQRT2 = np.sqrt(2)     # sqrt(2) with default precision np.float64


def hellinger1(p, q):
    return norm(np.sqrt(p) - np.sqrt(q)) / _SQRT2


def hellinger2(p, q):
    return euclidean(np.sqrt(p), np.sqrt(q)) / _SQRT2


def hellinger3(p, q):
    return np.sqrt(np.sum((np.sqrt(p) - np.sqrt(q)) ** 2)) / _SQRT2