I have fit a number of time series flux measurements of stars undergoing a physical process. Two of the parameters of these fits are important model parameters that tie back to the physics. I have samples from these fits and I intend to model the population distribution of these two model parameters following the method laid out in this paper. The likelihood function for this is given by:
My model is basically a 2d gaussian mixture model (these are the f’s in the sum, where alpha are the hyperparameters). I am fixing the mu and sigma for several of the 2d gaussians and only letting one float. I am sampling in the multiplicative factors and the mean vector and covariance matrix of the last 2d gaussian. The p_0 are the priors for the individual model fits, and K and N are the number of samples and number of stars, respectively. I’ve extracted the chains into a numpy array with shape (N,K,2) called X below.
I had trouble using Mvnormal, so I’ve settled on the following for my model and likelihood function:
@theano.compile.ops.as_op(itypes=[T.dvector,T.dscalar,T.dscalar,T.dscalar,T.dscalar,T.dscalar],otypes=[T.dscalar])
def loglike(phi_norm,tE_pbh,piE_pbh,sigma11_pbh,sigma22_pbh,sigma12_pbh):
g0 = scipy.stats.multivariate_normal.pdf(X, mean=mu0, cov=sigma0)
g1 = scipy.stats.multivariate_normal.pdf(X, mean= mu1, cov=sigma1)
g2 = scipy.stats.multivariate_normal.pdf(X, mean= mu2, cov=sigma2)
g3 = scipy.stats.multivariate_normal.pdf(X, mean= mu3, cov=sigma3)
mu = np.array([floating_mu0,floating_mu1])
cov = np.array([[sigma11,sigma12],[sigma12,sigma22]])
g4 = scipy.stats.multivariate_normal.pdf(tEpiE, mean=mu, cov=cov)
fa = phi_norm[0]*g0 + phi_norm[1]*g1 + phi_norm[2]*g2 + phi_norm[3]*g3 + phi_norm[4]g4
sums = 1/Knp.sum(fa/p0,axis=1)
return np.log10(np.prod(sums))
#pymc3 model
hierarchical_model = pm.Model()
with hierarchical_model:
#weight for mixture model, one for each gaussian
phi = pm.Uniform(‘phi’,lower=0,upper=N,shape=5)
# normalize them to the number of stars
phi_norm = pm.Deterministic(‘phi_norm’,N*phi/T.sum(phi))
# let the fifth gaussian float
floating_mu0 = pm.Normal(‘floating_mu0’,mu=np.log10(200),sd=2)
floating_mu1 = pm.Normal(‘floating_mu1’,mu=np.log10(0.01),sd=-1)
# sample each component of the covariance matrix (setting diagonal terms equal)
sigma11 = pm.Normal(‘sigma11’,mu=0.05,sd=0.1)
sigma22 = pm.Normal(‘sigma22’,mu=0.05,sd=0.1)
sigma12 = pm.Normal(‘sigma12’,mu=0.02,sd=0.05)
# send to likelihood function
pm.Potential(‘loglike’,loglike(phi_norm,floating_mu0, floating_mu1,sigma11,sigma22,sigma12))
trace = pm.sample(
tune=250,
draws=250,
chains=4,
step=pm.Metropolis(),
progressbar=True)
If I could figure out the MvNormal I could use NUTS instead of Metropolis by avoiding the @theano.compile.ops.as_op line above the likelihood function, and I could use help there. Even still, the code above runs into a chain failed error as below:
RemoteTraceback Traceback (most recent call last)
RemoteTraceback:
“”"
Traceback (most recent call last):
File “/anaconda3/lib/python3.7/site-packages/theano/compile/function_module.py”, line 903, in call
self.fn() if output_subset is None else
File “/anaconda3/lib/python3.7/site-packages/theano/gof/op.py”, line 892, in rval
r = p(n, [x[0] for x in i], o)
File “/anaconda3/lib/python3.7/site-packages/theano/compile/ops.py”, line 555, in perform
outs = self.__fn(*inputs)
File “”, line 10, in loglike
g4 = multivariate_normal.pdf(tEpiE, mean=mu, cov=cov)
File “/anaconda3/lib/python3.7/site-packages/scipy/stats/_multivariate.py”, line 521, in pdf
psd = _PSD(cov, allow_singular=allow_singular)
File “/anaconda3/lib/python3.7/site-packages/scipy/stats/_multivariate.py”, line 160, in init
raise ValueError(‘the input matrix must be positive semidefinite’)
ValueError: the input matrix must be positive semidefinite
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File “/anaconda3/lib/python3.7/site-packages/pymc3/parallel_sampling.py”, line 110, in run
self._start_loop()
File “/anaconda3/lib/python3.7/site-packages/pymc3/parallel_sampling.py”, line 160, in _start_loop
point, stats = self._compute_point()
File “/anaconda3/lib/python3.7/site-packages/pymc3/parallel_sampling.py”, line 191, in _compute_point
point, stats = self._step_method.step(self._point)
File “/anaconda3/lib/python3.7/site-packages/pymc3/step_methods/compound.py”, line 27, in step
point, state = method.step(point)
File “/anaconda3/lib/python3.7/site-packages/pymc3/step_methods/arraystep.py”, line 175, in step
apoint, stats = self.astep(self.bij.map(point))
File “/anaconda3/lib/python3.7/site-packages/pymc3/step_methods/metropolis.py”, line 162, in astep
accept = self.delta_logp(q, q0)
File “/anaconda3/lib/python3.7/site-packages/theano/compile/function_module.py”, line 917, in call
storage_map=getattr(self.fn, ‘storage_map’, None))
File “/anaconda3/lib/python3.7/site-packages/theano/gof/link.py”, line 325, in raise_with_op
reraise(exc_type, exc_value, exc_trace)
File “/anaconda3/lib/python3.7/site-packages/six.py”, line 692, in reraise
raise value.with_traceback(tb)
File “/anaconda3/lib/python3.7/site-packages/theano/compile/function_module.py”, line 903, in call
self.fn() if output_subset is None else
File “/anaconda3/lib/python3.7/site-packages/theano/gof/op.py”, line 892, in rval
r = p(n, [x[0] for x in i], o)
File “/anaconda3/lib/python3.7/site-packages/theano/compile/ops.py”, line 555, in perform
outs = self.__fn(*inputs)
File “”, line 10, in loglike
g4 = multivariate_normal.pdf(tEpiE, mean=mu, cov=cov)
File “/anaconda3/lib/python3.7/site-packages/scipy/stats/_multivariate.py”, line 521, in pdf
psd = _PSD(cov, allow_singular=allow_singular)
File “/anaconda3/lib/python3.7/site-packages/scipy/stats/_multivariate.py”, line 160, in init
raise ValueError(‘the input matrix must be positive semidefinite’)
ValueError: the input matrix must be positive semidefinite
Apply node that caused the error: FromFunctionOp{loglike}(phi_norm, tE_pbh_shared, piE_pbh_shared, sigma11_pbh_shared, sigma22_pbh_shared, InplaceDimShuffle{}.0)
Toposort index: 15
Inputs types: [TensorType(float64, vector), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar)]
Inputs shapes: [(5,), (), (), (), (), ()]
Inputs strides: [(8,), (), (), (), (), ()]
Inputs values: [array([712., 712., 712., 712., 712.]), array(2.30103), array(-2.), array(0.05), array(0.05), array(0.69606642)]
Outputs clients: [[MakeVector{dtype=‘float64’}(logp_phi_interval, __logp_tE_pbh, TensorConstant{-inf}, __logp_sigma11_pbh, __logp_sigma22_pbh, __logp_sigma12_pbh, loglike)]]
Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
File “/anaconda3/lib/python3.7/site-packages/ipykernel/zmqshell.py”, line 536, in run_cell
return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 2848, in run_cell
raw_cell, store_history, silent, shell_futures)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 2874, in _run_cell
return runner(coro)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/async_helpers.py”, line 67, in _pseudo_sync_runner
coro.send(None)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3049, in run_cell_async
interactivity=interactivity, compiler=compiler, result=result)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3214, in run_ast_nodes
if (yield from self.run_code(code, result)):
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3296, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File “”, line 16, in
pm.Potential(‘loglike’,loglike(phi_norm,tE_pbh,piE_pbh,sigma11_pbh,sigma22_pbh,sigma12_pbh))
HINT: Use the Theano flag ‘exception_verbosity=high’ for a debugprint and storage map footprint of this apply node.
“”"
The above exception was the direct cause of the following exception:
ValueError Traceback (most recent call last)
ValueError: the input matrix must be positive semidefinite
Apply node that caused the error: FromFunctionOp{loglike}(phi_norm, tE_pbh_shared, piE_pbh_shared, sigma11_pbh_shared, sigma22_pbh_shared, InplaceDimShuffle{}.0)
Toposort index: 15
Inputs types: [TensorType(float64, vector), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar), TensorType(float64, scalar)]
Inputs shapes: [(5,), (), (), (), (), ()]
Inputs strides: [(8,), (), (), (), (), ()]
Inputs values: [array([712., 712., 712., 712., 712.]), array(2.30103), array(-2.), array(0.05), array(0.05), array(0.69606642)]
Outputs clients: [[MakeVector{dtype=‘float64’}(logp_phi_interval, __logp_tE_pbh, TensorConstant{-inf}, __logp_sigma11_pbh, __logp_sigma22_pbh, __logp_sigma12_pbh, loglike)]]
Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
File “/anaconda3/lib/python3.7/site-packages/ipykernel/zmqshell.py”, line 536, in run_cell
return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 2848, in run_cell
raw_cell, store_history, silent, shell_futures)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 2874, in _run_cell
return runner(coro)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/async_helpers.py”, line 67, in _pseudo_sync_runner
coro.send(None)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3049, in run_cell_async
interactivity=interactivity, compiler=compiler, result=result)
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3214, in run_ast_nodes
if (yield from self.run_code(code, result)):
File “/anaconda3/lib/python3.7/site-packages/IPython/core/interactiveshell.py”, line 3296, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File “”, line 16, in
pm.Potential(‘loglike’,loglike(phi_norm,tE_pbh,piE_pbh,sigma11_pbh,sigma22_pbh,sigma12_pbh))
HINT: Use the Theano flag ‘exception_verbosity=high’ for a debugprint and storage map footprint of this apply node.
The above exception was the direct cause of the following exception:
RuntimeError Traceback (most recent call last)
in
6 chains=4,
7 step=pm.Metropolis(),
----> 8 progressbar=True)
9
/anaconda3/lib/python3.7/site-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, **kwargs)
435 _print_step_hierarchy(step)
436 try:
–> 437 trace = _mp_sample(**sample_args)
438 except pickle.PickleError:
439 _log.warning(“Could not pickle model, sampling singlethreaded.”)
/anaconda3/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, **kwargs)
967 try:
968 with sampler:
–> 969 for draw in sampler:
970 trace = traces[draw.chain - chain]
971 if (trace.supports_sampler_stats
/anaconda3/lib/python3.7/site-packages/pymc3/parallel_sampling.py in iter(self)
391
392 while self._active:
–> 393 draw = ProcessAdapter.recv_draw(self._active)
394 proc, is_last, draw, tuning, stats, warns = draw
395 if self._progress is not None:
/anaconda3/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
295 else:
296 error = RuntimeError(“Chain %s failed.” % proc.chain)
–> 297 raise error from old_error
298 elif msg[0] == “writing_done”:
299 proc._readable = True
RuntimeError: Chain 3 failed.