<commit_before>[Gibbs sampling](http://en.wikipedia.org/wiki/Gibbs_sampling) is a [Markov chain Monte Carlo](http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) method for sampling from a multivariate probability distribution. It has many applications in inference, [Bayesian networks](http://en.wikipedia.org/wiki/Bayesian_network) and machine learning. I am a bioinformatician, who needs to sample from joint distributions all the time. What motivated me to write this post is that currently I'm building a Bayesian network for my work, where I need to draw samples from the network. I'll talk about the use of Gibbs sampling in BayesNets in another post.

Let's start with a simple sampling problem. Suppose we have a single variable $latex X$ which takes two values:
<br>
\\( P(X =0)=0.5 \\) and \\( P(X =1)=0.5 \\)
<br>
This is an example of [binomial distribution](http://en.wikipedia.org/wiki/Binomial_distribution). How can get a sample from this distribution?
Simply, flip a coin. If it's head, \\( X=1 \\), else \\( X=0 \\). The following python code will generate samples from this distribution.

{% highlight python %}
import random

def sampler():
	r = random.random()
	if (r < 0.5):
		X = 0
	else:
		X = 1
	return X
{% endhighlight %}
<br>
What if you have a [multinomial distribution](http://en.wikipedia.org/wiki/Multinomial_distribution) ? Say you want to model the roll of a dice:
<br>
$$ latex P(X = i) = 1/6 ,  i \in 1,\ldots ,6$$
<br>
Divide the interval [0, 1] into 6 equal parts and select the value of \\( X \\) based on the interval in which a pseudo-random number falls.
```
def sampler():
    r = random.random()
    if (r < float(1)/6):
        X = 1
    elif (float(1)/6 <= r < float(2)/6):
        X = 2
    elif (float(2)/6 <= r < float(3)/6):
        X = 3
    elif (float(3)/6 <= r < float(4)/6):
        X = 4
    elif (float(5)/6 <= r < float(5)/6):
        X = 5
    else:
        X = 6
	return X
```
<br>
Suppose we have a multivariate distribution -
<br>
$$ P(X\_1,X\_2, \dots,X\_n) $$
<br>
If the variables are independent of each other (i.e [independence assumption](http://en.wikipedia.org/wiki/Independence_(probability_theory)), then we have
<br>
$$ P(X\_1,X\_2, \ldots,X\_n) = \prod\_{i=1}^{n} P(X\_i) $$
<br>
which will be easy to sample. Draw samples from each of the \\( P(X\_i) \\) individually and then multiply to get a sample of the joint distribution.
<br>
People who work with machine learning and data mining would recognize this assumption to be true for [Naive Bayes algorithm](http://en.wikipedia.org/wiki/Naive_Bayes_classifier).<commit_msg>gibbs sampling3<commit_after>[Gibbs sampling](http://en.wikipedia.org/wiki/Gibbs_sampling) is a [Markov chain Monte Carlo](http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) method for sampling from a multivariate probability distribution. It has many applications in inference, [Bayesian networks](http://en.wikipedia.org/wiki/Bayesian_network) and machine learning. I am a bioinformatician, who needs to sample from joint distributions all the time. What motivated me to write this post is that currently I'm building a Bayesian network for my work, where I need to draw samples from the network. I'll talk about the use of Gibbs sampling in BayesNets in another post.

Let's start with a simple sampling problem. Suppose we have a single variable $latex X$ which takes two values:
<br>
\\( P(X =0)=0.5 \\) and \\( P(X =1)=0.5 \\)
<br>
This is an example of [binomial distribution](http://en.wikipedia.org/wiki/Binomial_distribution). How can get a sample from this distribution?
Simply, flip a coin. If it's head, \\( X=1 \\), else \\( X=0 \\). The following python code will generate samples from this distribution.

{% highlight python %}
import random

def sampler():
	r = random.random()
	if (r < 0.5):
		X = 0
	else:
		X = 1
	return X
{% endhighlight %}
<br>
What if you have a [multinomial distribution](http://en.wikipedia.org/wiki/Multinomial_distribution) ? Say you want to model the roll of a dice:
<br>
$$ latex P(X = i) = 1/6 ,  i \in 1,\ldots ,6$$
<br>
Divide the interval [0, 1] into 6 equal parts and select the value of \\( X \\) based on the interval in which a pseudo-random number falls.
```python
def sampler():
    r = random.random()
    if (r < float(1)/6):
        X = 1
    elif (float(1)/6 <= r < float(2)/6):
        X = 2
    elif (float(2)/6 <= r < float(3)/6):
        X = 3
    elif (float(3)/6 <= r < float(4)/6):
        X = 4
    elif (float(5)/6 <= r < float(5)/6):
        X = 5
    else:
        X = 6
	return X
```
<br>
Suppose we have a multivariate distribution -
<br>
$$ P(X\_1,X\_2, \dots,X\_n) $$
<br>
If the variables are independent of each other (i.e [independence assumption](http://en.wikipedia.org/wiki/Independence_(probability_theory)), then we have
<br>
$$ P(X\_1,X\_2, \ldots,X\_n) = \prod\_{i=1}^{n} P(X\_i) $$
<br>
which will be easy to sample. Draw samples from each of the \\( P(X\_i) \\) individually and then multiply to get a sample of the joint distribution.
<br>
People who work with machine learning and data mining would recognize this assumption to be true for [Naive Bayes algorithm](http://en.wikipedia.org/wiki/Naive_Bayes_classifier).
{'fraction_non_alphanumeric': 0.1316193399846508, 'fraction_numerical': 0.018610897927858788, 'mean_word_length': 3.501727115716753, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2118754', 'n_tokens_mistral': 1746, 'n_tokens_neox': 1639, 'n_words': 659}