Bimodal Distribution – Recovering Parameters Using pymc3

bimodalnormal distributionpymc3python

I am trying to determine the parameters mu1, mu2, sigma1, sigma2, and w of a bimodal distribution using pymc3. x ~ w * Norm(u1, sigma1) + (1-w) * Norm(u1, sigma2)

I use the following code:

# Generate sample data
import numpy as np
from pylab import concatenate, normal

# First normal distribution parameters
mu1 = 1
sigma1 = 0.1

# Second normal distribution parameters
mu2 = 2
sigma2 = 0.2

w1 = 2/3 # Proportion of samples from first distribution
w2 = 1/3 # Proportion of samples from second distribution

n = 7500 # Total number of samples

n1 = int(n*w1) # Number of samples from first distribution
n2 = int(n*w2) # Number of samples from second distribution

# Generate n1 samples from the first normal distribution and n2 samples from the second normal distribution
y = concatenate((normal(mu1, sigma1, n1), normal(mu2, sigma2, n2)))


# Recover parameters mu1, sigma1, mu2, sigma2 and w from the sample data
import pymc3 as pm
from scipy.stats import norm
model = pm.Model()

def logp(mu, sigma, w, y):
  mu1 = mu[0]
  mu2 = mu[1]
  sigma1 = sigma[0]
  sigma2 = sigma[1]
  return np.log(w * norm.pdf(y, mu1, sigma1) + (1-w) * norm.pdf(y, mu2, sigma2)).sum()

with model:

  # Priors for unknown model parameters
  mu = pm.Normal("mu", mu=0, sigma=1, shape=2)
  sigma = pm.HalfNormal("sigma", sigma=1, shape=2)
  w = pm.distributions.continuous.Uniform("w")

  # Likelihood (sampling distribution) of observations
  likelihood=pm.DensityDist('Likelihood', logp, observed=dict(mu=mu, sigma=sigma, w=w, y=y))

start = pm.find_MAP(model = model)

The following is the error stack trace:

KeyError                                  Traceback (most recent call last)
/usr/local/lib/python3.7/dist-packages/theano/tensor/type.py in dtype_specs(self)
    279                 "complex64": (complex, "theano_complex64", "NPY_COMPLEX64"),
--> 280             }[self.dtype]
    281         except KeyError:

KeyError: 'object'

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
12 frames
<ipython-input-30-12971af8f9da> in <module>()
     25 
     26   # Likelihood (sampling distribution) of observations
---> 27   likelihood=pm.DensityDist('Likelihood', logp, observed=dict(mu=mu, sigma=sigma, w=w, y=y))
     28 
     29 start = pm.find_MAP(model = model)

/usr/local/lib/python3.7/dist-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
    120         else:
    121             dist = cls.dist(*args, **kwargs)
--> 122         return model.Var(name, dist, data, total_size, dims=dims)
    123 
    124     def __getnewargs__(self):

/usr/local/lib/python3.7/dist-packages/pymc3/model.py in Var(self, name, dist, data, total_size, dims)
   1165                     distribution=dist,
   1166                     total_size=total_size,
-> 1167                     model=self,
   1168                 )
   1169             self.observed_RVs.append(var)

/usr/local/lib/python3.7/dist-packages/pymc3/model.py in __init__(self, name, data, distribution, total_size, model)
   1871             datum.missing_values for datum in self.data.values() if datum.missing_values is not None
   1872         ]
-> 1873         self.logp_elemwiset = distribution.logp(**self.data)
   1874         # The logp might need scaling in minibatches.
   1875         # This is done in `Factor`.

<ipython-input-30-12971af8f9da> in logp(mu, sigma, w, y)
     15   sigma1 = sigma[0]
     16   sigma2 = sigma[1]
---> 17   return np.log(w * norm.pdf(y, mu1, sigma1) + (1-w) * norm.pdf(y, mu2, sigma2)).sum()
     18 
     19 with basic_model:

/usr/local/lib/python3.7/dist-packages/scipy/stats/_distn_infrastructure.py in pdf(self, x, *args, **kwds)
   1738         args = tuple(map(asarray, args))
   1739         dtyp = np.find_common_type([x.dtype, np.float64], [])
-> 1740         x = np.asarray((x - loc)/scale, dtype=dtyp)
   1741         cond0 = self._argcheck(*args) & (scale > 0)
   1742         cond1 = self._support_mask(x, *args) & (scale > 0)

/usr/local/lib/python3.7/dist-packages/theano/tensor/var.py in __truediv__(self, other)
    168 
    169     def __truediv__(self, other):
--> 170         return theano.tensor.basic.true_div(self, other)
    171 
    172     def __floordiv__(self, other):

/usr/local/lib/python3.7/dist-packages/theano/graph/op.py in __call__(self, *inputs, **kwargs)
    248         """
    249         return_list = kwargs.pop("return_list", False)
--> 250         node = self.make_node(*inputs, **kwargs)
    251 
    252         if config.compute_test_value != "off":

/usr/local/lib/python3.7/dist-packages/theano/tensor/elemwise.py in make_node(self, *inputs)
    497         using DimShuffle.
    498         """
--> 499         inputs = list(map(as_tensor_variable, inputs))
    500         out_dtypes, out_broadcastables, inputs = self.get_output_info(
    501             DimShuffle, *inputs

/usr/local/lib/python3.7/dist-packages/theano/tensor/basic.py in as_tensor_variable(x, name, ndim)
    205         )
    206 
--> 207     return constant(x, name=name, ndim=ndim)
    208 
    209 

/usr/local/lib/python3.7/dist-packages/theano/tensor/basic.py in constant(x, name, ndim, dtype)
    253         assert x_.ndim == ndim
    254 
--> 255     ttype = TensorType(dtype=x_.dtype, broadcastable=[s == 1 for s in x_.shape])
    256 
    257     try:

/usr/local/lib/python3.7/dist-packages/theano/tensor/type.py in __init__(self, dtype, broadcastable, name, sparse_grad)
     52         # True or False
     53         self.broadcastable = tuple(bool(b) for b in broadcastable)
---> 54         self.dtype_specs()  # error checking is done there
     55         self.name = name
     56         self.numpy_dtype = np.dtype(self.dtype)

/usr/local/lib/python3.7/dist-packages/theano/tensor/type.py in dtype_specs(self)
    281         except KeyError:
    282             raise TypeError(
--> 283                 f"Unsupported dtype for {self.__class__.__name__}: {self.dtype}"
    284             )
    285 

TypeError: Unsupported dtype for TensorType: object

Best Answer

The following code from the pymc3 docs would be a better way to run a normal mixture model:

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt

print(f"Running on PyMC3 v{pm.__version__}")

np.random.seed(12345)  # set random seed for reproducibility

k = 3
ndata = 500
spread = 5
centers = np.array([-spread, 0, spread])

# simulate data from mixture distribution
v = np.random.randint(0, k, ndata)
data = centers[v] + np.random.randn(ndata)

plt.hist(data);

# setup model
model = pm.Model()
with model:
    # cluster sizes
    p = pm.Dirichlet("p", a=np.array([1.0, 1.0, 1.0]), shape=k)
    # ensure all clusters have some points
    p_min_potential = pm.Potential("p_min_potential", tt.switch(tt.min(p) < 0.1, -np.inf, 0))
    # cluster centers
    means = pm.Normal("means", mu=[0, 0, 0], sigma=15, shape=k)
    # break symmetry
    order_means_potential = pm.Potential(
        "order_means_potential",
        tt.switch(means[1] - means[0] < 0, -np.inf, 0)
        + tt.switch(means[2] - means[1] < 0, -np.inf, 0),
    )
    # measurement error
    sd = pm.Uniform("sd", lower=0, upper=20)
    # latent cluster of each observation
    category = pm.Categorical("category", p=p, shape=ndata)
    # likelihood for each observed value
    points = pm.Normal("obs", mu=means[category], sigma=sd, observed=data)


# fit model
with model:
    step1 = pm.Metropolis(vars=[p, sd, means])
    step2 = pm.ElemwiseCategorical(vars=[category], values=[0, 1, 2])
    tr = pm.sample(10000, step=[step1, step2], tune=5000)