MissingInputError: Undeclared input - Bimodal distribution

I am trying to determine what is wrong with the following code that aims to recover the parameters of a bimodal normal distribution. The code first generates some sample data with known parameters mu1, sigma1, mu2, sigma2, w1 and w2 then the code using pymc3 tries to get back the parameters:

from pylab import concatenate, normal
import pymc3 as pm
import numpy as np

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 = 500 # Total number of samples

n1 = int(n*w1) # Number of samples from first distribution
n2 = n-n1 # 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(
    'blah', logp, observed=dict(mu=mu, sigma=sigma, w=w, y=y))

with model:
step = pm.Metropolis(vars=[mu, sigma, w])
tr = pm.sample(1000, step=step, tune=500)

Running the code fails at the last line of code with the following stack trace:

/usr/local/lib/python3.7/dist-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
637 if idata_kwargs:
638 ikwargs.update(idata_kwargs)
→ 639 idata = arviz.from_pymc3(trace, **ikwargs)
641 if compute_convergence_checks:

/usr/local/lib/python3.7/dist-packages/arviz/data/io_pymc3.py in from_pymc3(trace, prior, posterior_predictive, log_likelihood, coords, dims, model, save_warmup, density_dist_obs)
570 model=model,
571 save_warmup=save_warmup,
→ 572 density_dist_obs=density_dist_obs,
573 ).to_inference_data()

/usr/local/lib/python3.7/dist-packages/arviz/data/io_pymc3.py in init(self, trace, prior, posterior_predictive, log_likelihood, predictions, coords, dims, model, save_warmup, density_dist_obs)
170 self.density_dist_obs = density_dist_obs
→ 171 self.observations, self.multi_observations = self.find_observations()
173 def find_observations(self) → Tuple[Optional[Dict[str, Var]], Optional[Dict[str, Var]]]:

/usr/local/lib/python3.7/dist-packages/arviz/data/io_pymc3.py in find_observations(self)
182 elif hasattr(obs, “data”) and self.density_dist_obs:
183 for key, val in obs.data.items():
→ 184 multi_observations[key] = val.eval() if hasattr(val, “eval”) else val
185 return observations, multi_observations

/usr/local/lib/python3.7/dist-packages/theano/graph/basic.py in eval(self, inputs_to_values)
552 inputs = tuple(sorted(inputs_to_values.keys(), key=id))
553 if inputs not in self._fn_cache:
→ 554 self._fn_cache[inputs] = theano.function(inputs, self)
555 args = [inputs_to_values[param] for param in inputs]

/usr/local/lib/python3.7/dist-packages/theano/compile/function/init.py in function(inputs, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input)
348 on_unused_input=on_unused_input,
349 profile=profile,
→ 350 output_keys=output_keys,
351 )
352 return fn

/usr/local/lib/python3.7/dist-packages/theano/compile/function/pfunc.py in pfunc(params, outputs, mode, updates, givens, no_default_updates, accept_inplace, name, rebuild_strict, allow_input_downcast, profile, on_unused_input, output_keys)
530 profile=profile,
531 on_unused_input=on_unused_input,
→ 532 output_keys=output_keys,
533 )

/usr/local/lib/python3.7/dist-packages/theano/compile/function/types.py in orig_function(inputs, outputs, mode, accept_inplace, name, profile, on_unused_input, output_keys)
1976 on_unused_input=on_unused_input,
1977 output_keys=output_keys,
→ 1978 name=name,
1979 )
1980 with config.change_flags(compute_test_value=“off”):

/usr/local/lib/python3.7/dist-packages/theano/compile/function/types.py in init(self, inputs, outputs, mode, accept_inplace, function_builder, profile, on_unused_input, fgraph, output_keys, name)
1582 # make the fgraph (copies the graph, creates NEW INPUT AND
→ 1584 fgraph, additional_outputs = std_fgraph(inputs, outputs, accept_inplace)
1585 fgraph.profile = profile
1586 else:

/usr/local/lib/python3.7/dist-packages/theano/compile/function/types.py in std_fgraph(input_specs, output_specs, accept_inplace)
186 orig_outputs = [spec.variable for spec in output_specs] + updates
→ 188 fgraph = FunctionGraph(orig_inputs, orig_outputs, update_mapping=update_mapping)
190 for node in fgraph.apply_nodes:

/usr/local/lib/python3.7/dist-packages/theano/graph/fg.py in init(self, inputs, outputs, features, clone, update_mapping)
161 for output in outputs:
→ 162 self.import_var(output, reason=“init”)
163 for i, output in enumerate(outputs):
164 self.clients[output].append((“output”, i))

/usr/local/lib/python3.7/dist-packages/theano/graph/fg.py in import_var(self, var, reason)
340 "Computation graph contains a NaN. " + var.type.why_null
341 )
→ 342 raise MissingInputError(“Undeclared input”, variable=var)
343 self.setup_var(var)
344 self.variables.add(var)

MissingInputError: Undeclared input

Your observed argument isn’t pointing to (only) data and your logp() function is not (directly) related to your model parameters. It’s not entirely clear to me why you are using pm.DensityDist here. I suspect you want pm.Potential. I would suggest this fully-worked example of performing inference on a mixture of Gaussians for some pointers.

I have seen the fully worked example. However, what I am trying to do is learn how to create my own custom distribution by providing a logpdf. In my understanding, observed can either be used to send the observed data or a dictionary containing the data and named parameters to be used by logpdf.

In that case, I would suggest looking at this example in which a mixture model is put together “by hand” with a custom likelihood using pm.Potential. Jun Peng’s talk on pymc3 likelihoods may also be useful. He goes through some pretty custom likelihoods, but he mostly uses the plumbing provided by pym3.

I found this lecture note from Standford titled (Gaussian mixture models). Equation (1) bottom of page 10 shows the log likelihood that I am using. Top of page 11 states: “Direct maximization of L(θ; z) is quite difficult numerically, because of the sum of terms inside the logarithm.” It goes on to explain how to use a latent variable to solve the problem quite similarly to what the fully worked example you provided. In fact, I got my code running but it was not converging to any meaningful solution, probably falling into some local maxima.
Thanks for your time helping me!