I am trying to determine the parameters mu1, mu2, sigma1, sigma2, and w of a bimodal distribution using pymc3. 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
Updated the code as per accepted answer. Crash no longer happens but I do not get a convergence towards the solution:
import pymc3 as pm
import numpy as np
from pylab import concatenate, normal
# Generate sample data
# 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)))
def normpdf(y, mu, sigma):
return (1 / sigma) * pm.math.exp(pow(y - mu, 2) / pow(sigma, 2))
def logp(mu, sigma, w, y):
mu1 = mu[0]
mu2 = mu[1]
sigma1 = sigma[0]
sigma2 = sigma[1]
w1 = w[0]
w2 = w[1]
return pm.math.sum(np.log(w1 * normpdf(y, mu1, sigma1) + w2 * normpdf(y, mu2, sigma2)))
# Recover parameters mu1, sigma1, mu2, sigma2 and w from the sample data
model = pm.Model()
with model:
# Priors for unknown model parameters
mu = pm.Normal("mu", mu=np.array([0.5, 1.5]), sigma=1, shape=2)
sigma = pm.HalfNormal("sigma", sigma=1, shape=2)
w = pm.Dirichlet("w", a=np.array([1.1, 0.9]), shape=2)
# 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)
print(start)
Output:
{'mu': array([-2.12059455, 1.95271142]),
'sigma': array([0.01268331, 0.37437335]),
'sigma_log__': array([-4.36746845, -0.98250172]),
'w': array([0.70625152, 0.29374848]),
'w_stickbreaking__': array([0.43862378])}
Using pdb I stepped through and found that the error was on this line:
return np.log(w * norm.pdf(y, mu1, sigma1) + (1-w) * norm.pdf(y, mu2, sigma2)).sum()
because scipy.stats.norm doesn't know how to deal with pymc/Theano objects. While I'm sure there's an easy way to get a normal density from pymc3 objects, it's also straightforward to define a function to calculate it ourselves:
def normpdf(y, mu, sigma):
return pm.math.sum((1 / sigma) * pm.math.exp(pow(y - mu, 2) / pow(sigma, 2)))
def logp(mu, sigma, w, y):
mu1 = mu[0]
mu2 = mu[1]
sigma1 = sigma[0]
sigma2 = sigma[1]
return np.log(w * normpdf(y, mu1, sigma1) + (1-w) * normpdf(y, mu2, sigma2)).sum()
# Recover parameters mu1, sigma1, mu2, sigma2 and w from the sample data
model = pm.Model()
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")
w = pm.Dirichlet("p", a=np.array([1.0, 1.0]), shape=2)
# 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)
This should work; I've not inspected the output but it at least compiles and runs without error.