Gibbs Sampler of Bivariate Normal

Let's say that we want to generate samples from a bivariate normal with the following parameter.
\[\mu=\left[ {\begin{array}{*{20}{c}}{{\mu _1}}\\{{\mu _2}}\\\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}1\\2\\\end{array}} \right]\]
\[\Sigma=\left[ {\begin{array}{*{20}{c}}1 & \rho\\\rho & 1 \\\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}1 & {0.9} \\{0.9} & 1 \\\end{array}} \right]\]
The conditional distribution $f(x_1|x_2)$ is univariate normal $N(\mu_1+\rho(x_2-\mu_2), 1-\rho^2)$ and $f(x_2|x_1)$ is $N(\mu_2+\rho(x_1-\mu_1), 1-\rho^2)$.

The Matlab code for the Gibbs sampling is

n=5000;
x=zeros(n,2);
rho=0.9;
mu=[1;2];
sig=sqrt(1-rho^2);
%initial point
x(1,:)=[10,10];
for i=2:n   
    x(i,1)=normrnd(mu(1)+rho*(x(i-1,2)-mu(2)),sig);
    x(i,2)=normrnd(mu(2)+rho*(x(i,1)-mu(1)),sig);
end

and the result after 1,000 samples burn-in is

plot(x(1001:end,1),x(1001:end,2),'.');
axis square;

To check the result, calculate the sample mean and sample covariance matrix of x

>> mean(x)
ans =
    0.9942    1.9960
>> corrcoef(x)
ans =
    1.0000    0.8943
    0.8943    1.0000

This problem is, in fact, easily solved with Matlab's mvnrnd().

covm=[1 rho;rho 1];
X = mvnrnd(mu,covm,4000);
figure;
plot(X(:,1),X(:,2),'.');
axis square;

Conjugate Prior

Before I said anything about the title, recall Mr. Baye's rule
\[p(\theta |y) = \frac{{p(y|\theta )}p(\theta)}{{p(y)}}\]
The denominator $p(y)$ is called marginal , and is treated as a constant, therefore
\[p(\theta |y) \propto p(y|\theta )p(\theta )\]
The $p(\theta)$ is the prior. The $p(y|\theta)$, which should be familiar by now, is the ... likelihood, and $p(\theta|y)$ is the posterior.
To get the posterior, we must calculate the product of the likelihood and the prior.

Now, we might collect some data and calculate the $\theta$ for the likelihood. What if the likelihood is of a super complex model with ugly parameters? You will end up with the even uglier posterior.

Instead of doing that, we do the trick. Look for a way to choose $p(\theta)$ so that $p(\theta|y)$ has the same function form as $p(\theta)$. That is the idea of Conjugate Prior.

For example (the not so ugly likelihood), If your likelihood is Poisson
\[p(y|\theta ) = \frac{{{e^{ - \theta }}{\theta ^y}}}{{y!}}\]
The nice conjugate prior is Gamma, or $\theta \sim G(\alpha ,\beta )$

\[p(\theta ) = \frac{{{\theta ^{\alpha - 1}}{e^{ - \theta /\beta }}}}{{\Gamma (\alpha ){\beta ^\alpha }}}\]
Therefore,
\[p(\theta |y) \propto \left( {{e^{ - \theta }}{\theta ^y}} \right)({\theta ^{\alpha - 1)}}{e^{ - \theta /\beta }})\]
\[ = {\theta ^{y + \alpha - 1}}{e^{ - \theta (1 + 1/\beta )}}\]
The formula for the posterior is essentially the same as the prior with some rescaling, in fact
\[\theta |y \sim G(y + \alpha ,\frac{1}{{1 + 1/\beta }})\]

Fitting Data with MATLAB's "dfittool"

If we have data and a distribution in our mind, "dfittool" can help us find the distribution's parameters and visualize the result.

Let's generate normal samples (mu=5, sigma=2) in MATLAB and fit it with a normal distribution  using "dfittool".

 >> x = normrnd(5,2,100,1);

That's very good fit, isn't it? Looking at the results, it finds the mu to be 4.98 and sigma 2.05, which can be calculated from mean and standard deviation of x.

>> mu = mean(x)
mu =
    4.9822

>> sig = sqrt(var(x))
sig =
    2.0518

It also says "Log Likelihood: -213.267". Where does it come from?

Recall that the pdf of normal distribution is
\[f(X|(\mu,{\sigma ^2})) = \frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{{{{(x - \mu)}^2}}}{{2{\sigma ^2}}}}}\]
and, therefore, the likelihood function is
\[\varphi (\mu,{\sigma ^2}) = \prod\limits_{i = 1}^n {\frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{{{{({x_i} - \mu)}^2}}}{{2{\sigma ^2}}}}}} \]
and the log-likelihood function is
\[\log \varphi (\mu,{\sigma ^2}) = \sum\limits_{i = 1}^n {(\log \frac{1}{{\sqrt {2\pi } }} - \log \sigma - \frac{{{{({x_i} - \mu)}^2}}}{{2{\sigma ^2}}})} \]
\[ = n\log \frac{1}{{\sqrt {2\pi } }} - n\log \sigma - \frac{1}{{2{\sigma ^2}}}\sum\limits_{i = 1}^n {{{({x_i} - \mu )}^2}} \]
Substitute mu and sigma to the last equation, using the MATLAB command 
>> n*log(1/sqrt(2*pi))-n*log(sig)- (1/(2*sig^2))*sum((x-mu).^2)
ans =
 -213.2672

Maximum Likelihood Estimation

In statistics, there are two schools of thoughts, the frequentist and the Bayesian (more about them in later blogs). One of the most popular tools for the frequentist is Maximum Likelihood Estimations (MLE).

 \[\hat \theta = \mathop {\arg \max }\limits_{all\theta } L(\theta |{x_1},...x{}_n)\]

The procedure goes like this.

1. Collect the data. For example, height, weight, etc.
2. Make an assumption about the distribution of the data, i.e. Normal Distribution (the shape of the likelihood $f({x_1},...,{x_n}|\mu ,\sigma)$), but with unknown mean and variance.  The mean and variance are parameters we want to estimate.
3. Let's assume they are many ways to get the mean and variance, i.e. different values of mu and sigma. Which mean value gives highest likelihood? also Which variance value give hightest likelihood?
4. It turns out that, for this case, the answers are sample mean $\hat\mu = {\textstyle{{\sum {{x_i}} } \over n}}$, and sample variance $\hat\sigma = \sqrt {\frac{{\sum {({x_i} - \hat\mu} {)^2}}}{n}} $ that we are so familiar with.

See the proof here.

Another example: A coin was tossed 80 times, yields 49 heads, and 31 tails but we do not know which coin. There are three coins in the basket with probability of giving head 1/3, 1/2 and 2/3. Which coin was used for this activity?

Naive answer: The proportion of the heads is  49/80 = 0.6125. Three coins have the head prob = [0.3333 0.5 0.6666]. The 0.6125 is closer to 0.6666 than any other values. Therefore, we think coin #3 must be the one.

Better answer: The distribution of the data is binomial. With different values of p, we have:
>> binopdf(49,80,1/3)
ans =
  2.0789e-007

>> binopdf(49,80,1/2)
ans =
    0.0118

>> binopdf(49,80,2/3)
ans =
    0.0545

Coin #3 has maximum likelihood, therefore it must be coin #3.

The next question, What is the MLE of p?

>> simplify(diff(p^49 * (1-p)^31))
ans =
-p^48*(80*p - 49)*(p - 1)^30

To make the above equation equal to zero (and therefore maximize the likelihood), the answer for p is 0, 1 and 49/80. Thus MLE for p is 49/80 (same as the number of heads in the data).

Gibbs Sampler

Gibbs sampler is a special case of Metropolis-Hasting algorithm. The differences are that first, we always accepts a candidate point and second, we must know the full conditional distributions.

Let's assume that we have a joint distribution $f(x,y_1,...y_d)$, we would like to know more about the marginal density $f(x)$. So, we must integrate over all other variables.
\[f(x) = \int {...\int {f(x,{y_1},...,{y_d})d{y_1}...d{y_d}} } \]
As with other MCMC methods, Gibbs sampler generates $X_1,X_2,...,X_m$ from $f(x)$ and then use the sample to estimate the desired characteristic of $f(x)$.

For bivariate case, to generate samplers (variates) from $f(x,y)$ we alternately sampling  $x_i \sim f(x|y_{i - 1})$ and $y_i \sim f(y|x_i)$ until we got the number of samples we want.

Example. Consider the following joint distribution
\[f(x,y) \propto \left({\begin{array}{*{20}{c}}n\\x\\\end{array}} \right) y^{x+\alpha-1}(1 - y)^{n-x+\beta-1}\]
where $ x = 0,1,...,n $ and $0 \le y \le 1$.
This is the joint distribution of a binomial distribution
\[f(x|n,y)= \left({\begin{array}{*{20}{c}}n\\x\\\end{array}} \right) y^x(1 - y)^{n-x}\]
and a beta distribution
\[f(y|\alpha ,\beta ) = \frac{1}{B(\alpha ,\beta )}y^{\alpha-1}(1-y)^{\beta-1}\]
(confused? Try multiplying these two distributions together, you got the idea)

In short, \[X|y \sim Bino(n,y)\] and \[Y|x \sim Beta(x+\alpha,n-x+\beta)\]

The MATLAB code to generate samples using the Gibbs Sampler is:

% starting points;
x(1)=binornd(n,0.5,1,1);
y(1)=betarnd(x(1)+a,n-x(1)+b,1,1);
for i = 2:k
    x(i)=binornd(n,y(i-1),1,1);
    y(i)=betarnd(x(i)+a,x-x(i)+b,1,1);
end

The estimation from the Gibbs sampler is very close to the true marginal (Casella & George 1992):

\[f(x) = \left( {\begin{array}{*{20}{c}}n\\x\\\end{array}} \right)\frac{{\Gamma (\alpha+\beta)\Gamma (x+\alpha)\Gamma (n-x+\beta )}}{{\Gamma (\alpha )\Gamma (\beta )\Gamma (\alpha+\beta+n)}}\]

The code for plotting the left figure is just the histogram of $X$

[c,xi]=hist(x(m+1:end),17);
bar(xi,c/sum(c),1);

and on the right figure, we compute the true marginal density

for x = 0:16
    fhat(x+1)=nchoosek(n,x)*gamma(a+b)*gamma(x+a)*gamma(n-x+b)/...
        (gamma(a)*gamma(b)*gamma(a+b+n));
end
bar(xi,fhat,1);

Metropolis Sampler

The original Metropolis algorithm (1953) is different from the Metropolis-Hasting algorithm (1970), where the proposal distribution is symmetric (i.e. normal distribution). Thus,
\[q(Y|X)=q(X|Y)\]
and yield
\[\alpha(X_t,Y)=min \left\{1,\frac{\pi(Y)}{\pi(X_t)}\right\}\]
when the proposal distribution is such that $q(X|Y)=q(Y|X)=q(|X-Y|)$, then it is also called "random walk Metropolis"

Metropolis-Hasting Algorithm

One efficient way to get samples from complex distribution is by using Markov Chain. Metropolis-Hasting Algorithm is a famous one for this purpose.

The candidate point for the next sample is accepted as the next state of the chain with the probability given by

\[\alpha (X_t,Y) = \min \left\{ {1,\frac{\pi (Y)q(X_t|Y)}{\pi (X_t)q(Y|X_t)} \right\}\]
If the point $Y$ is not accepted, the chain does not move and  $X_{t+1}=X_t$

$q(.|X_t)$ is a proposal distribution. It can be multivariate normal with mean $X_t$ and fixed covariance matrix.

For example, generate random variables from a standard Cauchy distribution (target distribution) given by

\[f(x) = \frac{1}{\pi} \left[\frac{1}{(1 + {x^2})}\right]; \;\:\:-\infty < x < \infty \]
or
\[f(x) \propto \frac{1}{{(1 + {x^2})}}\]

We can use a normal distribution for the proposal distribution.

The MATLAB code for Metropolis-Hasting algorithm is as follow.

n = 10000;

sig = 2;

x = zeros(1,n);

x(1)=normrnd(0,sig);

cauchy = @(x) 1./(1+x.^2);

norm = @(x,mu,sig) 1/sig * exp(-0.5*((x-mu)/sig).^2);

for i=2:n  

    y = normrnd(x(i-1),sig);

    alpha = min([1, cauchy(y)*norm(x(i-1),y,sig)/...

                   (cauchy(x(i-1))*norm(y,x(i-1),sig))]);

    if rand() < alpha

        x(i)=y;

    else

        x(i)=x(i-1);

    end

end

%discard the first 500 points

cv = x(501:end);

% Plot the distribution and compare to the true cauchy distribution

[n,x]=hist(cv,50);

width=x(2)-x(1);

bar(x,n/(sum(n)*width),1); hold on;

lx = linspace(min(x),max(x),100);

plot(x,cauchy(x)/pi,'r','Linewidth',2); hold off;

axis tight;