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)
640
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()
574
/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)
169
170 self.density_dist_obs = density_dist_obs
→ 171 self.observations, self.multi_observations = self.find_observations()
172
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
186
/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]
556
/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 )
534
/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
1583 # OUTPUT VARIABLES)
→ 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
187
→ 188 fgraph = FunctionGraph(orig_inputs, orig_outputs, update_mapping=update_mapping)
189
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)
160
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