1
$\begingroup$

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
$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Community Bot
    Sep 1, 2021 at 8:58
  • $\begingroup$ @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. $\endgroup$
    – Tarik
    Sep 1, 2021 at 10:44
  • $\begingroup$ @Community This question stackoverflow.com/questions/35990467/… uses GaussianMixture to solve this. I would like to solve the same problem using pymc3 library. $\endgroup$
    – Tarik
    Sep 1, 2021 at 10:45
  • $\begingroup$ You probably should't cross-post $\endgroup$ Sep 1, 2021 at 18:00

2 Answers 2

1
$\begingroup$

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)
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Tarik
    Sep 1, 2021 at 17:08
  • $\begingroup$ 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. $\endgroup$
    – Tarik
    Sep 2, 2021 at 7:59
0
$\begingroup$

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"]);
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.