# Recovering Bimodal distribution parameters using pymc3

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:

52         # True or False
---> 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

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Sep 1, 2021 at 8:58
• @Community I have sample data that I want to model as a bimodal distribution (statisticshowto.com/what-is-a-bimodal-distribution). I want to get the parameters mu1, mu2, sigma1, sigma2 and w1 & w2 of the binomial distribution best fitting the sample data. The first part of the provided code generates sample data. The second part attempts to determine the distribution parameters. Sep 1, 2021 at 10:44
• @Community This question stackoverflow.com/questions/35990467/… uses GaussianMixture to solve this. I would like to solve the same problem using pymc3 library. Sep 1, 2021 at 10:45
• You probably should't cross-post Sep 1, 2021 at 18:00

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

# 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)

• Thanks for the time you spent answering my question. I will carefully study your code (I am a noob in statistics). That said, would it be possible to have a look at what I have done and explain what's wrong with it? Sep 1, 2021 at 17:08
• OK, I have modified the code to have separate sd values and it worked like a charm. It's a much better way to get things done indeed. Sep 2, 2021 at 7:59

The following is the modified code as per alan ocallaghan advice:

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

# 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

ndata = 5000 # Total number of samples

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

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

k = 2

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

plt.hist(data);

# setup model
model = pm.Model()
with model:
# cluster sizes
p = pm.Dirichlet("p", a=np.array([1.0] * k), 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=list(range(k)), sigma=15, shape=k)

# break symmetry
order_means_potential = pm.Potential("order_means_potential", tt.switch(means[1] - means[0] < 0, -np.inf, 0))

# measurement error
sd = pm.Uniform("sd", lower=0, upper=20, shape=2)

# 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[category], observed=data)

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

pm.traceplot(tr, var_names=["p", "sd", "means"]);