PyMC Tutorial #3: Latent Dirichlet Allocation (LDA) using PyMC

Before you read this post, we suggest you to read our previous post regarding Naïve Bayes NB topic model since the code presented in this post is just the modification from the previous post. Some terms have been explained in the previous post.

Previously, we have explained that Naïve Bayes (NB) model assumes that each document in the collection is drawn from a single topic. Unlike Naïve Bayes, Latent Dirichlet Allocation (LDA) assumes that a single document is a mixture of several topics [1][2]. Moreover, we will use smoothed version of LDA, which is described in its original paper authored by Blei et al. [1]. The smoothed version of LDA was proposed to tackle the sparsity problem in our collection. A document might have words that do not appear in the other documents. As a result, this situation will set zero probability to such document. Hence, it is impossible to generate the document. Smoothed-LDA overcomes this kind of problem [1].

The “generative aspect” of LDA model

We will just modify a little bit the NB model from the previous post !

Suppose, \(D\) is the number of documents in our collection. \(N_d\) is the number of words in the \(d\)-th document. \(K\) is the number of predefined-topics. In this case, the number of topics is not automatically inferred. Instead, we will manually set the value of \(K\) based on our intuition.

For each topic \(k\), we draw its word distribution, which is denoted as \(\phi_k\). Since we are going to model the word distribution using multinomial/categorical distribution, we then use the Dirichlet distribution (its conjugate prior) to model its parameter.

\[\phi_k \sim Dir(\beta), 1 \leq k \leq K\]

After that, for each document \(d\), we draw a topic distribution, which is denoted as \(\theta_d\)

\[\theta_d \sim Dir(\alpha), 1 \leq d \leq D\]

See ! the preceding notation differs from NB model ! since a document might contain more than one topics (with their corresponding proportions) !

The hyperparameters \(\alpha\) and \(\beta\) are assumed to be fixed in our case.

Then, for each word \(n\) in document \(d\), we draw a topic for that word. This topic is denoted as \(z_{d,n}\)

\[z_{d,n} \sim Categorical(\theta_d), 1 \leq d \leq D, 1 \leq n \leq N_d\]

In NB model, we have seen that a single topic is associated with a document. On the other hand, LDA associates a single topic with a single word. Hence, drawing topic must be performed in the level of word, not document.

Finally, after we know the topic of a word (\(z_{d,n}\)), we then draw the physical word itself from the word distribution associated with its selected topic. Each word is denoted as \(w_{d,n}\).

\[w_{d, n} \sim Categorical(\phi_k),  1 \leq d \leq D, 1 \leq n \leq N_d\]

The following figure is the plate notation for LDA model (smoothed):

Now, it is time to build the code using PyMC. We just need to modify the previous code :)

-code1-

 import numpy as np  
 import pymc as pc  
   
   
 def wordDict(collection):  
  word_id  = {}  
  idCounter = 0  
  for d in collection:  
    for w in d:  
      if (w not in word_id):  
        word_id[w] = idCounter  
        idCounter+=1  
  return word_id  
   
 def toNpArray(word_id, collection):  
  ds = []  
  for d in collection:  
    ws = []  
    for w in d:  
      ws.append(word_id.get(w,0))  
    ds.append(ws)  
  return np.array(ds)  
   
 ###################################################  
   
 #doc1, doc2, ..., doc7  
 docs = [["sepak","bola","sepak","bola","bola","bola","sepak"],  
         ["uang","ekonomi","uang","uang","uang","ekonomi","ekonomi"],  
         ["sepak","bola","sepak","bola","sepak","sepak"],  
         ["ekonomi","ekonomi","uang","uang"],  
         ["sepak","uang","ekonomi"],  
         ["komputer","komputer","teknologi","teknologi","komputer","teknologi"],  
         ["teknologi","komputer","teknologi"]]  
   
 word_dict = wordDict(docs)  
 collection = toNpArray(word_dict,docs)  
   
 #number of topics  
 K = 3  
   
 #number of words (vocab)  
 V = len(word_dict)  
   
 #number of documents  
 D = len(collection)  
   
 #array([1, 1, 1, ..., 1]) K times  
 alpha = np.ones(K)  
   
 #array([1, 1, 1, ..., 1]) V times  
 beta = np.ones(V)  
   
 #array containing the information about doc length in our collection
 Nd = [len(doc) for doc in collection]  
   
   
 ######################## LDA model ##################################  
   
 #topic distribution per-document  
 theta = pc.Container([pc.CompletedDirichlet("theta_%s" % i,   
                                             pc.Dirichlet("ptheta_%s"%i, theta=alpha))  
                      for i in range(D)])  
   
 #word distribution per-topic  
 phi = pc.Container([pc.CompletedDirichlet("phi_%s" % j,   
                                           pc.Dirichlet("pphi_%s" % j, theta=beta))  
                     for j in range(K)])  
   
   
 #Please note that this is the tricky part :)  
 z = pc.Container([pc.Categorical("z_%i" % d,  
                                  p = theta[d],  
                                  size = Nd[d],  
                                  value = np.random.randint(K, size=Nd[d]))   
                   for d in range(D)])  
   
 #word generated from phi, given a topic z  
 w = pc.Container([pc.Categorical("w_%i_%i" % (d,i),  
                                  p = pc.Lambda("phi_z_%i_%i" % (d,i),  
                                                lambda z=z[d][i], phi=phi : phi[z]),
                                  value=collection[d][i],  
                                  observed=True)  
                   for d in range(D) for i in range(Nd[d])])  
   
 ####################################################################  
   
 model = pc.Model([theta, phi, z, w])  
 mcmc = pc.MCMC(model)  
 mcmc.sample(iter=5000, burn=1000)  
   
   
 #show the topic assignment for each word, using the last trace  
 for d in range(D):  
    print(mcmc.trace('z_%i'%d)[3999])  
   

In our computer, the output of the preceding code is as follows:

-code2-

 [2 2 2 2 2 2 2]  
 [1 1 1 0 0 1 1]  
 [2 2 2 2 2 2]  
 [1 0 0 0]  
 [1 0 1]  
 [1 2 0 1 1 1]  
 [1 1 1]  
   

What does it mean ? Let’s compare the preceding results with our document collection:

-code3-

 docs = [["sepak","bola","sepak","bola","bola","bola","sepak"],  
         ["uang","ekonomi","uang","uang","uang","ekonomi","ekonomi"],  
         ["sepak","bola","sepak","bola","sepak","sepak"],  
         ["ekonomi","ekonomi","uang","uang"],  
         ["sepak","uang","ekonomi"],  
         ["komputer","komputer","teknologi","teknologi","komputer","teknologi"],  
         ["teknologi","komputer","teknologi"]]  
   

Remember that we set \(K = 3\). This value may come from our intuition, or our prior knowledge and experience.

We can see that the word “sepak” in doc #1 is associated with topic 2. The word “ekonomi” in doc #2 is associated with topic 1. The word “uang” in doc #5 is associated with topic 0. Using the results, we we can say that topic 2 is talking about sport, while topic 1 is somehow talking about economics. And, we cannot be really certain what topic 0 tells us.

We can also read the results as follows: doc #1 consists of only topic 2. Doc #2 is a mixture of topic 1 and topic 0. Doc #6 is generated from 3 topics.

We really cannot judge whether or not this LDA model works well here. Like we said before, it will be interesting if we can try this on a real word dataset to see the behavior :)

Main References:

[1] Latent Dirichlet Allocation, David M. Blei, Andrew Ng, Michael I. Jordan. Journal of Machine Learning Research, 2003.

[2] Integrating Out Multinomial Parameters in Latent Dirichlet Allocation and Naïve Bayes for Collapsed Gibbs Sampling, Bob Carpenter, LingPipe Inc, 2010.

Alfan Farizki Wicaksono

(firstname [at] cs [dot] ui [dot] ac [dot] id)

Fakultas Ilmu Komputer, UI

Ditulis di Depok, 27 Juli 2015

PyMC Tutorial #2: Estimating the Parameters of A Naïve Bayes (NB) Model

Before you read this post, we suggest you to read our previous post regarding PyMC Introduction, several parameter estimation techniques, as well as Bernoulli parameter estimation.

NB model can be considered as one of the topic models (like Latent Dirichlet Allocation (LDA)), in which each word in a document is generated based on topic-specific multinomial/categorical distribution [1]. In addition to that, NB assumes that each document is drawn from a single topic (NOT mixture of topics) [1].

The “generative aspect” of NB model

Suppose, \(D\) is the number of documents in our collection. \(N_d\) is the number of words in the \(d\)-th document. \(K\) is the number of predefined-topics. In this case, the number of topics is not automatically inferred. Instead, we will manually set the value of \(K\) based on our intuition.

For each topic \(k\), we draw its word distribution, which is denoted as \(\phi_k\). Since we are going to model the word distribution using multinomial/categorical distribution, we then use the Dirichlet distribution (its conjugate prior) to model its parameter.

\[\phi_k \sim Dir(\beta),   1 \leq k \leq K\]

After that, we draw a global (collection) topic distribution, which is denoted as \(\theta\).

\[\theta \sim Dir(\alpha)\]

The hyperparameters \(\alpha\) and \(\beta\) are assumed to be fixed in our case.

Then, for each document \(d\) we draw a topic describing that document. This topic is denoted as \(z_m\).

\[z_m \sim Categorical(\theta),   1 \leq d \leq D\]

Here, Instead of modeling \(z_m\) as a multinomial random variable, we model it using Categorical distribution. We do not have any particular reason, except that later it will be easy for us to build our code this way.

Finally, after we know the topic of a document, we then draw words from the word distribution associated with its selected topic. Each word is denoted as \(w_{d,n}\).

\[w_{d, n} \sim Categorical(\phi_k),  1 \leq d \leq D, 1 \leq n \leq N_d\]

The following figure describes the plate notation of NB model

Estimating marginal posterior distribution of \(\theta\)

Our goal is to estimate \(\theta\), that is, we want to automatically infer the topic for each document, which is currently hidden. All we know is just the data, that is, all documents with their words. In other words, our task is just like a clustering task ! we aim at clustering documents that have the same topic.

First, we import the needed library for this problem.

-code1-

 import numpy as np  
 import pymc as pc  
 import random as rn  

Suppose, for our toy-implementation, our document collection is represented as follows.

-code2-

 #doc1, doc2, ..., doc7  
   
 docs = [["sepak","bola","sepak","bola","bola","bola","sepak"],  
         ["uang","ekonomi","uang","uang","uang","ekonomi","ekonomi"],  
         ["sepak","bola","sepak","bola","sepak","sepak"],  
         ["ekonomi","ekonomi","uang","uang"],  
         ["sepak","uang","ekonomi"],  
         ["komputer","komputer","teknologi","teknologi","komputer","teknologi"],  
         ["teknologi","komputer","teknologi"]]  

Then, we need to convert such collection into a special form such that we will be easy to process the collection later. In this case, we will create a set of words, and associate each word with a unique integer (as an id). Instead using array of words, we will use array of id’s.

-code3-

 def wordDict(collection):  
  word_id  = {}  
  idCounter = 0  
  for d in collection:  
    for w in d:  
      if (w not in word_id):  
         word_id[w] = idCounter  
         idCounter+=1  
  return word_id  
   
 def toNpArray(word_id, collection):  
  ds = []  
  for d in collection:  
     ws = []  
     for w in d:  
        ws.append(word_id.get(w,0))  
     ds.append(ws)  
  return np.array(ds)  

After that, we need to set several useful variables. Suppose, we set the number of topic as 3.

-code4-

 word_dict = wordDict(docs)  
 collection = toNpArray(word_dict,docs)  
   
 #number of topics  
 K = 3  
   
 #number of words (vocab)  
 V = len(word_dict)  
   
 #number of documents  
 D = len(collection)  
   
 #array([1, 1, 1, ..., 1]) K times  
 alpha = np.ones(K)  
   
 #array([1, 1, 1, ..., 1]) V times  
 beta = np.ones(V)  
   
 #array containing the information about doc length in our collection  
 Nd = [len(doc) for doc in collection]  

Now, let us see how easy for us if we want to create our own topic model using PyMC :) the following code will explain how to represent NB model using PyMC.

-code5-

 #naive bayes model  
   
 #word distribution for each topic  
 phi = pc.Container([pc.CompletedDirichlet("phi_%i" % k,   
                      pc.Dirichlet("pphi_%i" % k, theta=beta))  
           for k in range(K)])  
   
 #topic distribution in the collection  
 theta = pc.CompletedDirichlet("theta",   
                               pc.Dirichlet("ptheta", theta=alpha))  
   
 #for each document, draw a topic z_m  
 z = pc.Container([pc.Categorical("z_%i" % d,  
                                  p = theta,  
                                  value = rn.randint(0,K-1))  
                  for d in range(D)])  
   
 #for each document, draw words, based on topic z_m  
 w = pc.Container([pc.Categorical("w_%i_%i" % (d,i),  
                                  p = pc.Lambda("phi_z_%i_%i" % (d,i),  
                                                lambda z=z[d], phi=phi : phi[z]),  
                                  value=collection[d][i],  
                                  observed=True)  
                  for d in range(D) for i in range(Nd[d])])  
   
 model = pc.Model([theta, phi, z, w])  

“Container” class is used when we need to wrap all the identical random variables (such as \(\phi_k\)'s) in a single reference. “Dirichlet” is a PyMC class for Dirichlet distribution. Unfortunately, it only generates \(k-1\) values since the \(k\)-th value can be easily determined by the fact that the sum of all probability values should be \(1\). Hence, “CompletedDirichlet” class will add the \(k\)-th value to the last part. “Lambda” is a class in PyMC if we want to represent a deterministic random variable using lambda expression. You can read more about all these stuffs via PyMC homepage.

Next, let’s perform MCMC to estimate the parameter of NB model !

-code6-

 mcmc = pc.MCMC(model)  
 mcmc.sample(iter=5000, burn=1000)  

In the preceding code, we use 5000 iterations (to generate the sample as well as improve the quality of marginal posterior distribution of the parameter), and assuming that burn-in period is achieved after we perform 1000 iterations. The generated traces will not consider the results before burn-in period.

Finally, let us see the estimation results ! Remember once again, that previously we used \(K = 3\), which means that we are really sure that there are 3 topics in our collection, but we do not know which document has which topic !

-code7-

 #show z_i after sampling (random walks) has ended  
 #number of traces = iter – burn = 5000 – 1000 = 4000  
   
 for d in range(D):  
    print(mcmc.trace('z_%i'%d)[3999])  

In our computer, the results are as follows:

-code8-

 1  
 2  
 1  
 2  
 1  
 0  
 0  

Let us compare the results with our collection once again:

-code9-

 docs = [["sepak","bola","sepak","bola","bola","bola","sepak"],                  #doc1
         ["uang","ekonomi","uang","uang","uang","ekonomi","ekonomi"],            #doc2
         ["sepak","bola","sepak","bola","sepak","sepak"],                        #doc3
         ["ekonomi","ekonomi","uang","uang"],                                    #doc4
         ["sepak","uang","ekonomi"],                                             #doc5
         ["komputer","komputer","teknologi","teknologi","komputer","teknologi"], #doc6
         ["teknologi","komputer","teknologi"]]                                   #doc7

We can see that doc #1, #3, and #5 are assigned the same topic (topic 1).  Doc #2 and #4 are of the same topic as well (topic 2). Lastly, Doc #6 and #7 have the same topic (topic 0). We can see that this algorithm works pretty well on automatically clustering the documents within the same topic. However, this is just a toy-implementation. We need to be really careful about this result.

It will be interesting if we can try this program in a real-world document collection :)

Main References:

[1] Integrating Out Multinomial Parameters in Latent Dirichlet Allocation and Naïve Bayes for Collapsed Gibbs Sampling, Bob Carpenter, LingPipe Inc, 2010.

Alfan Farizki Wicaksono

(firstname [at] cs [dot] ui [dot] ac [dot] id)

Fakultas Ilmu Komputer, UI

Ditulis di Depok, 27 Juli 2015

PyMC Tutorial #1: Bayesian Parameter Estimation for Bernoulli Distribution

Suppose we have a Coin which consists of two sides, namely Head (H) and Tail (T). All of you might know that we can model a toss of a Coin using Bernoulli distribution, which takes the value of \(1\) (if H appears) with probability \(\theta\) and \(0\) (if T appears) with probability \(1 - \theta\). In this case, \(\theta\) is also called as the parameter of a Bernoulli distribution since knowing the value of \(\theta\) is sufficient for determining \(P(H)\) and \(P(T)\). For a fair Coin, \(\theta\) is set to \(0.5\), which means that we have equal degree of belief for both sides.

This time, we aim at estimating the parameter \(\theta\) of a particular Coin. To do that, first, we need to collect the data sample, which serves as our evidence, from an experiment. Second, we use that data to estimate the parameter \(\theta\). Suppose, to collect the data, we toss the Coin 10 times and record the outcomes. We get a sequence of \(\{H, H, T, H, ..., T\}\) which consists of 10 elements, in which each element represents the outcome a single coin tossing. By assuming that the previous data sample is independent and identically distributed (often referred to as i.i.d), we then perform statistical computation to determine the estimate of \(\theta\).

There are two broad categories of estimating the parameter of a known probability distribution. The first one is so called Maximum Likelihood Estimation (MLE) and the second one is Bayesian parameter estimation. We will examine both methods briefly in this post. In the end, we will focus on Bayesian parameter estimation and show the usage of PyMC (Python library for MCMC framework) to estimate the parameter of a Bernoulli distribution.


Maximum Likelihood Estimation (MLE)

Please do not be afraid when you hear the name of this method ! Eventhough the name of this method is somewhat “long-and-complicated”, but the opposite situation actually happens. MLE often involves basic counting of events on our data. As an example, MLE estimates the paramater θ of the Coin using the following, “surprisingly simple”, statistic

\[\hat{\theta} = \frac{\# Heads}{\# Heads + \# Tails}\]

Because of that, people usually refer MLE as a “Frequentist approach”. In general, MLE aims at seeking a set of parameters which maximizes the likelihood of seeing our data.

\[\hat{\theta} = \substack{argmax \\ \theta} P(x_1, x_2, ..., x_n|\theta)\]

Now, let us try to implement MLE for estimating the parameter of a Bernoulli distribution (using Python programming language). We simulate the experiment of tossing a Coin N times using a list of integer values, in which 1 and 0 represents Head and Tail, respectively. Each value is generated randomly from a Bernoulli distribution.

\[P(H) = P(1) = \theta\]
\[P(T) = P(0) = 1 - \theta\]

We use Bernoulli-like distribution provided by Scipy library. So, we need to import this library as the first step.

-code1-

 from scipy.stats import bernoulli  

Next, we generate a sample data using the following code.

-code2-

 sampleSize = 20  
 theta = 0.2  
   
 def generateSample(t, s):  
  return bernoulli.rvs(t, size=s)  
   
 data = generateSample(theta, sampleSize)  

The preceding code will assign “data” with the following value

-code3-

 array([1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0])  

We can see that assigning theta to 0.2 makes the number of 0’s much more than the number of 1’s, which means that Tail has higher probability to occur compared to Head.

Now, we pretend that we do not know the parameter \(\theta\) and we only know the data. Given that data, we are going to estimate the value of \(\theta\), which is unknown to us. We use MLE, which means that we need to implement the aforementioned statistic

-code4-

 def thetaMLE(data):  
   count = 0  
   for i in data:  
      count+=i  
   return count/float(len(data))  

Now, let see several estimates when we use different sample size.

-code5-

 def showSeveralEstimates(sampleSizes):  
  for size in sampleSizes:  
    estimate = thetaMLE(generateSample(0.2, size))  
    print("using sample with size %i : theta = %f" % (size,estimate))  
   
 showSeveralEstimates([10,100,1000,2000,5000,10000])  

The preceding code yields the following results (the results may differ each time you run this program since it involves random sampling):

-code6-

 using sample with size 10 : theta = 0.3  
 using sample with size 100 : theta = 0.23  
 using sample with size 1000 : theta = 0.194  
 using sample with size 2000 : theta = 0.1965  
 using sample with size 5000 : theta = 0.1982  
 using sample with size 10000 : theta = 0.2006  

Look ! we can see that as the size of data increase, the estimate is getting closer to the real value of \(\theta\), that is \(0.2\). This concludes that if you want to obtain better estimate, you need to increase your data.


Bayesian Parameter Estimation

Although MLE is often easy to prepare as well as to compute, it has several limitations. One of them is that MLE can not leverage prior information or knowledge regarding the parameter itself. For example, based on our experience, we are really certain that a Coin is fair. Unfortunately, when we try to estimate the parameter using MLE, we cannot incorporate such knowledge to our computation. On the other hand, Bayesian Parameter Estimation takes into account prior knowledge regarding the parameter, which makes Bayesian Parameter Estimation provides more realistic and accurate estimation [1][2]. Sometimes we have prior belief about something before the observance of the data or evidence. But, once we finally see the evidence or data about that, we may change our belief [1][2].

Instead of directly estimating \(P(data|parameter)\), Bayesian Parameter Estimation estimates \(P(parameter|data)\).  Here, prior information about parameter \(\theta\) is encoded as a probability distribution \(P(\theta)\), which means that we consider \(\theta\) as a value of a random variable. When we quantify uncertainty about \(\theta\), it will be easy for us to encode our prior belief. After we observe our data, we then change our prior belief towards \(\theta\), into our posterior belief, denoted as \(P(\theta|X)\).

\[P(\theta|x_1, x_2, ..., x_n) \propto P(\theta) P(x_1, x_2, ..., x_n|\theta) \]

In the preceding formula, \(P(\theta)\) is the prior distribution of \(\theta\). \(P(X|\theta)\) is the likelihood of our observed data. The likelihood represents how likely that we will see the observed data when we know already the parameter \(\theta\). \(P(\theta |X)\) is the posterior distribution that represents the belief about \(\theta\) after taking both the data and prior knowledge into account (after we see our data).

We usually use the expected value to give the best estimate of \(\theta\). In other words, given the data \(X\), the estimate of \(\theta\) is obtained by calculating \(E[\theta |X]\). \(E[\theta |X]\) is usually called as the Bayes estimator.

\[\hat{\theta} = E[\theta |x_1, x_2, ..., x_n] \]

Hierarchical Bayesian Model

The prior distribution \(P(\theta)\) may be estimated using the so called hyperprior distributions. This kind of model is known as Hierarchical Bayesian Model. Furthermore, we can also estimate the hyperprior distribution itself, using hyper-hyperprior distribution, and so on. The reason behind using hyperprior distribution is that, instead we use directly the distribution \(\theta\), which may be available (from previous experiment), why don’t we let the “present data tell us about \(\theta\) by themselves” [2].

Let us see the previous example, in which we try to estimate the parameter of Bernoulli parameter \(\theta\), given the data collected by conducting several tosses of a Coin \(\{H, T, H, H, H, T,  ..., T\}.\) Suppose \(x_i\) represents the value of a single Coin toss.

\[x_i \sim Ber(\theta)\]

Now, we can model the parameter \(\theta\) using Beta distribution. In other words, \(\theta\) is a random variable that follows Beta distribution with parameter \(\alpha\) and \(\beta\). \(\alpha\) and \(\beta\) are called hyper-parameters. We use Beta distribution since it is the prior conjugate of a Bernoulli distribution. We will not elaborate more on the notion of conjugacy in this post. However, there are several mathematical reasons behind the use of conjugate prior distributions. One of them is that conjugate prior distribution will make our computation easier.

\[\theta \sim Beta(\alpha, \beta)\]

The posterior distribution of \(\theta\) can be then denoted as follows

\[P(\theta |x_1, ..., x_n, \alpha, \beta) \propto P(\theta |\alpha, \beta) P(x_1, ..., x_n |\theta, \alpha, \beta)\]

We can also represent the preceding model using the well-known plate notation as follows

Where \(N\) represents the number of tosses that we perform (the size of sample data). We get back to our main goal: estimating \(\theta\) (the posterior distribution of \(\theta\)) using Bayesian parameter estimation. We have just learnt that the estimation task involves computing the expectation value of \(\theta\) (\(E[\theta |X]\)), which means that we might need to perform a number of integrations. Unfortunately, in some cases, performing integrations will not be feasible, or at least it will be difficult to achieve a specified accuracy. Thus, we need to think of any approximation methods to back up our plan.

There are many types of numerical approximations for Bayesian parameter estimation. One of them (the most common) is Markov Chain Monte Carlo (MCMC). MCMC estimates the posterior distribution of \(\theta\) by performing a number of iterations or sampling. In each iteration, we improve the quality of our target distribution by leveraging the sampled data, and hoping that it will eventually arrive at the “true” posterior distribution of \(\theta\).

PyMC: A Python Library for MCMC Framework

Now, we are ready to play with the programming problem. Python has a library that provides MCMC framework for our problem. This library is called PyMC. You can go directly to its official website, if you want to know more about it.

First, let’s import several libraries that we need, including PyMC and pymc.Matplot for drawing histogram.

-code7-

 import pymc as pc  
 import pymc.Matplot as pt  
 import numpy as np  
 from scipy.stats import bernoulli  

Next, we need to create our model.

-code8-

 def model(data):  
   
   theta_prior = pc.Beta('theta_prior', alpha=1.0, beta=1.0)  
    
   coin = pc.Bernoulli('coin', p=theta_prior, value=data, observed=True)  
   
   mod = pc.Model([theta_prior, coin])  
   
   return mod  

In the preceding code, we represent \(\theta\) as “theta_prior”, which follows Beta distribution with parameter \(\alpha\) and \(\beta\). Here, we set both \(\alpha\) and \(\beta\) with 1.0. “coin” represents a sequence of coin tosses (NOT a single toss), in which each toss follows Bernoulli distribution (This corresponds to \(X\) in the preceding plate notation). We set “observed=True” since this is our observed data. “p=theta_prior” means that the parameter of “coin” is “theta_prior”. Here, our goal is to estimate the expected value of “theta_prior”, which is unknown. MCMC will perform several iterations to generate the sample from “theta_prior”, in which each iteration will improve the quality of the sample. Finally, we wrap all of our random variables using class Model.

Like the previous one, we need a modul that can generate our toy-sample:

-code9-

 def generateSample(t, s):  
   return bernoulli.rvs(t, size=s)  

Suppose, we have already generated a sample, and pretend that we do not know the parameter of the distribution where it comes from. We then use the generated sample to estimate \(\theta\).

-code10-

 def mcmcTraces(data):  
    
   mod = model(data)  
   mc = pc.MCMC(mod)  
    
   mc.sample(iter=5000, burn=1000)  
   return mc.trace('theta_prior')[:]  

The preceding procedure/function will produce traces, or MCMC samples generated by a number of interations. Based on that code, MCMC will iterate 5000 times. “burn” specifies the minimum iteration that we need before we are sure that we have achieved the “true” posterior distribution of \(\theta\). The function yields the traces of MCMC (except the sample generated before the burn-in period).

Now, let’s perform the MCMC run on our model, and plot the posterior distribution of \(\theta\) on a histogram.

-code11-

 sample = generateSample(0.7, 100)  
    
 trs = mcmcTraces(sample)  
 pt.histogram(trs, "theta prior; size=100", datarange=(0.2,0.9))  

Suppose the data is generated from a Bernoulli distribution with parameter \(\theta = 0.7\) (size = 100). If we draw the traces of \(\theta\) using Histogram, we will get the following figure.

We can see that the distribution is centered in the area 0.65 – 0.80. We are most likely happy with this result (since the prediction somehow close to 0.70), yet the variance is still very high. Now, let see what will happened when we increase the size of our observable data !

The following histogram was generated when we set size to 500:

The following histogram was generated when we set size to 5000:

See ! when we increase the size of our data, then the variance of the distribution gets lower. Thus, we are more confident about our prediction !

If we need the estimate value of \(\theta\), we can use the expected value (mean) of that distribution. We can use numpy library to get the mean of our sample.

-code12-

 #estimated theta  
   
 est_theta = np.mean(trs)  
 print(est_theta)  

Main References:
[1] Building Probabilistic Graphical Models with Python, Kiran R. Karkera, PACKT Publishing 2014
[2] Bayesian Inference, Byron Hall (STATISTICAT, LLC)

Alfan Farizki Wicaksono

(firstname [at] cs [dot] ui [dot] ac [dot] id)

Fakultas Ilmu Komputer, UI

Ditulis di Tambun, Bekasi, 26 Juli 2015