Possibly naive question, but is there a way to do structure learning of Bayesian networks with PyMC3? One possibility might be to define a score based on the overall model logp and then run an optimizer to maximize / minimize that score. The logp, I imagine, can be directly obtained from pymc3?
Thoughts? Other approaches?