How to do hidden variable learning in Bayesian Network with Python?
While this model could be implementable in the libpgm library (it seems to have quite rigid interface tailored to several special models, though), it will not allow you to reproduce the results of this paper. In the paper authors actually do something way simpler than Bayesian Inference – they only perform maximum-aposteriori (MAP) inference, which is way easier, but does tell you how confident the model is.

You don't have to code optimization algorithms for MAP inference by yourself, these days there are powerful and general packages that do almost everything for you. You can use automatic differentiation features of TensorFlow / PyTorch – they will compute the gradients and run optimization algorithm for you. You can also use numpy.autograd to get gradients for free, and then run any optimizer you like (check out scipy.optimize). The autograd + scipy.optimize pair is particularly neat if you wanna use more advanced optimization methods, e.g. 2nd order ones.

All these methods will give you values (point estimates) of $w, \theta, z$ that you can use to make inferences. However, if you want to have more than point estimates, and would like to do find all possible values (along with their probabilities) and maybe Bayesian Model Averaging, you'll need more powerful libraries enabling probabilistic programming. There are Edward (built upon TensorFlow), Pyro (built upon PyTorch), Stan and others.

If you feel overwhelmed by the choices, I recommend going with PyTorch.