# Naive Bayes model with PyMC3

**URL:** <https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314>\
**Category:** Questions\
**Created:** [December 4, 2018, 6:42am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314 "2018-12-04T06:42:00Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 4, 2018, 6:42am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/1 "2018-12-04T06:42:00Z")

</div>

Hello,

I’m new to PyMC3 and for getting started I tried to implement a Naive Bayes (NB) model as described in [this blog post](http://alfan-farizki.blogspot.com/2015/07/pymc-tutorial-2-estimating-parameters.html) but for version 3.5 of PyMC. Here’s my initial attempt:

```python
import numpy as np
import pymc3 as pm

K = 3 # number of topics
V = 4 # size of vocabulary

# 13 documents with 5 words each
data = np.array([[0, 0, 0, 0, 0], 
                 [0, 0, 0, 0, 1], 
                 [0, 0, 0, 1, 1], 
                 [0, 0, 0, 0, 1], 
                 [0, 0, 0, 0, 1], 
                 [0, 0, 0, 0, 1], 
                 [0, 0, 0, 0, 0], 
                 [1, 1, 1, 2, 2], 
                 [1, 1, 1, 1, 2], 
                 [1, 1, 1, 2, 1], 
                 [1, 1, 1, 2, 2], 
                 [2, 2, 2, 3, 3], 
                 [2, 2, 2, 3, 3]])

# number of documents
D = data.shape[0]

# Hyper-parameters
alpha = np.ones(K)
beta = np.ones(V)

with pm.Model() as model:
    # Global topic distribution
    theta = pm.Dirichlet("theta", a=alpha)
    
    # Word distributions for K topics
    phi = pm.Dirichlet("phi", a=beta, shape=(K, V))
    
    # Topic of documents
    z = pm.Categorical("z", p=theta, shape=D)
    
    # Words in documents
    for i in range(D):
        pm.Categorical(f"w_{i}", 
                       p=phi[z[i]], 
                       observed=data[i])

```

The results (not shown here) look reasonable but I’m curious if this implementation is correct at all and if there’s a better way to implement that with PyMC3, especially the `for` loop at the end of the example. I’m aware of more efficient alternatives for estimating the parameters of a NB model but I thought this is a good example for getting started with PyMC3.

Thanks in advance for any help!

Best Regards,  
Martin

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 4, 2018, 10:37am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/2 "2018-12-04T10:37:59Z")

</div>

I think your problem is quite similar to [latent dirichlet allocation](https://docs.pymc.io/notebooks/lda-advi-aevb.html) - you can have a look at this doc.

Otherwise, a general suggestion on mixture model like this kind is _not_ to explicitly model the latent categories (i.e., avoid doing `z = pm.Categorical("z", p=theta, shape=D)`), but rewrite it as a marginalized mixture model instead.

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 4, 2018, 11:31am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/3 "2018-12-04T11:31:26Z")

</div>

Thanks for your response. I’ve seen the document you mentioned but wanted to have a much simpler example running. The implementation was actually inspired by the [LDA example you posted earlier](https://discourse.pymc.io/t/lda-implementation-with-pymc3/1515/4).

Regarding

> [@junpenglao](#):
>
> … rewrite it as a marginalized mixture model instead

do you mean something comparable to [Marginalized Gaussian Mixture Model](https://docs.pymc.io/notebooks/marginalized_gaussian_mixture_model.html) but using the general [`Mixture`](https://docs.pymc.io/api/distributions/mixture.html#pymc3.distributions.mixture.Mixture) class instead?

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 5, 2018, 7:16am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/4 "2018-12-05T07:16:16Z")

</div>

I tried using the `Mixture` class following the examples in the [source code](https://github.com/pymc-devs/pymc3/blob/b40029969eb00dc88a791fb13b8076a8bfb02604/pymc3/distributions/mixture.py#L44) but got stuck applying it to my example here. Can you please show me how to re-write my example to a marginalized mixture model? I tried several things but had problems matching shapes of variables and observed data.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 5, 2018, 8:06am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/5 "2018-12-05T08:06:09Z")

</div>

Yes you do need to be careful to make the shape work… Maybe rewriting the observed likelihood function as a potential or Densitydist is a bit easier, since you can reuse the loglike function from the doc and the LDA post.

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 6, 2018, 5:23am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/6 "2018-12-06T05:23:49Z")

</div>

To verify my understanding I implemented a marginalized multivariate GMM. First, I generated some data :

```python
from sklearn.datasets import make_blobs

K = 5
n_samples = 150

X, y = make_blobs(n_samples=n_samples, centers=K)
plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.viridis, alpha=.25);

```

![fig-1](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/1/189ac6d575305b241ed167531d1c345cdd8384a3.png)

and then fit the following model to data `X` (shape=(150, 2)):

```python
cov = np.array([[1, 0], 
                [0, 1]])

with pm.Model() as model:
    w = pm.Dirichlet('w', np.ones(K))

    components = []

    for i in range(K):
        mu = pm.Normal('mu%i'%i, 0, 10, shape=2)
        components.append(pm.MvNormal.dist(mu=mu, cov=cov))

    xobs = pm.Mixture('x_obs', w, components, observed=X)

```

This works as expected, no problems here.

The marginalized version of the NB model I initially posted has the same structure (except for variable types):

```python
# Definitions of K, V, D, data as in my initial post
components = []

with pm.Model() as model:
    theta = pm.Dirichlet("theta", a=alpha)

    for i in range(K):
        phi = pm.Dirichlet(f"phi_{i}", a=beta)
        components.append(pm.Categorical.dist(p=phi))

    obs = pm.Mixture('obs', w=theta, comp_dists=components, observed=data)        

```

but in this case I get

```auto
ValueError: Input dimension mis-match. (input[0].shape[2] = 3, input[1].shape[2] = 5)

```

which I do not understand. I can get a working model definition when I define a `Mixture` _for each document_:

```python
with pm.Model() as model:
    theta = pm.Dirichlet("theta", a=alpha)

    for i in range(K):
        phi = pm.Dirichlet(f"phi_{i}", a=beta)
        components.append(pm.Categorical.dist(p=phi))

    for i in range(D):
        pm.Mixture(f'obs_{i}', w=theta, comp_dists=components, observed=data[i])

```

but in this case the sampling results don’t make sense. For example, inferred `theta` looks like

![fig-2](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/6/692683024b46e72b74396b4bb53043b483d57f0a.png)

I really would like to understand why

- `obs = pm.Mixture('obs', w=theta, comp_dists=components, observed=data)` leads to a dimension mis-match error
- fitting with the last model definition gives results for `theta` that do not make sense given the data. (With the non-marginalized NB model version from my initial post I get reasonable values for `theta`).

Please help me to answer these questions!

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 6, 2018, 9:13am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/7 "2018-12-06T09:13:22Z")

</div>

I will try using `DensityDist` later as you suggested but still would like to know why the previous approach using `Categorical` doesn’t work. I mean the example problem is so simple and I really would like to understand why this can’t be solved with built-in distributions/functions without resorting to `DensityDist`.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 6, 2018, 10:39am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/8 "2018-12-06T10:39:52Z")

</div>

I will try to have a look soon - a bit busy this week 😅

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 6, 2018, 11:02am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/9 "2018-12-06T11:02:36Z")

</div>

No rush, thanks a lot, appreciate it!

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 8, 2018, 7:19am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/10 "2018-12-08T07:19:46Z")

</div>

So here is how I would approach it:  
First, I would rewrite the explicit model, so that the observed is a flatten array - it is easier to handle shape wise:

```python
# flatten and index data
data_flatten = np.reshape(data, data.shape[0]*data.shape[1])
data_index = np.repeat(np.arange(data.shape[0]), data.shape[1])

with pm.Model() as model:
    # Global topic distribution
    theta = pm.Dirichlet("theta", a=alpha)
    
    # Word distributions for K topics
    phi = pm.Dirichlet("phi", a=beta, shape=(K, V))
    
    # Topic of documents
    z = pm.Categorical("z", p=theta, shape=D)
    
    # Words in documents
    p = phi[z][data_index]
    w = pm.Categorical("w", p=p, observed=data_flatten)
    trace = pm.sample(1000, tune=1000)

```

Now, to turn this into a marginalized model, we need to _not_ do indexing of the `p` in `p = phi[z]`, but keep the dimensions of the latent topics. Computationally, we want the logp to be evaluated for each latent topic (not just for the “true” latent label):

```python
with pm.Model() as model_marg:
    # Word distributions for K topics
    phi = pm.Dirichlet("phi", a=beta, shape=(K, V))
    
    # Topic of documents
    z = pm.Dirichlet("z", a=alpha/2., shape=(D, K))
    
    # Global topic distribution
    theta = pm.Deterministic("theta", z.mean(axis=0))
    
    # Words in documents
    comp_dists = pm.Categorical.dist(phi)
    w = pm.Mixture("w",
                   w=z[data_index, :],
                   comp_dists=comp_dists,
                   observed=data_flatten)
    trace_marg = pm.sample(1000, tune=1000)

```

you can check the mixture component shape to make sure that is the case:

```python
w.distribution.w.tag.test_value.shape
w.distribution._comp_logp(data_flatten).tag.test_value.shape

```

Of course, the two model is not identical, and there are inference problem like multi-modality for both of them… Mixture model is difficult unfortunately.

I put everything in a notebook here: [Planet\_Sakaar\_Data\_Science/discourse\_2314.ipynb at main · junpenglao/Planet\_Sakaar\_Data\_Science · GitHub](https://github.com/junpenglao/Planet_Sakaar_Data_Science/blob/master/PyMC_QnA/discourse_2314.ipynb)

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 9, 2018, 8:02am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/11 "2018-12-09T08:02:55Z")

</div>

Thanks a lot @junpenglao for taking time writing this up. It helped me a lot!

Regarding inference problems, I found that when using sparse Dirichlet priors i.e.

```python
alpha = np.ones(K) * 0.5
beta = np.ones(V) * 0.5

```

results significantly improve. For example, for the first model, inferred `theta` matches much better the actual global topic distribution:

 ![fig-1](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/2/26caebe49491027629b9b55c69edbf1e95cc9e83.png)

I couldn’t see that improvement for the marginalized mixture model though.

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [December 9, 2018, 8:24am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/12 "2018-12-09T08:24:54Z")

</div>

More informative prior almost always helps 😉

Using more sparse Dirichlet priors also help in the marginal case I think - you can see the mode switching more clearly.

---

<div class="post-metadata">

**Author:** ![krasserm](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/krasserm/32/1292_2.png) [@krasserm](https://discourse.pymc.io/u/krasserm)\
**Post date:** [December 9, 2018, 8:53am UTC](https://discourse.pymc.io/t/naive-bayes-model-with-pymc3/2314/13 "2018-12-09T08:53:19Z")

</div>

I used sparse Dirichlet priors here as this is [assumed for LDA](https://en.wikipedia.org/wiki/Latent_Dirichlet_allocation#Topics) (which is similar to the example here).