Maybe my description of the question is misleading. I do not want to do observed deterministic inference, but try to write codes to build customized likelihood from a Normal and a Bernoulli variable by
logp_ = w1 * ln( LL(feat1) ) + w2 * ln( LL(feat2) )
Where feat1 ~ Bernoulli distribution and feat2 ~ Normal distribution
Some one told me about using Poisson zeros trick (Bernoulli ones trick). But I want to know how to implement it in pymc3 (and want to know how these tricks work).
I am new to python and pymc3, it would be good to have an example that provide detailed annotations to help me learning how to do it and how to expand it to more than two random variables.
Really appreciate if you can help on this.
The following are my complete codes.
When running them, I got errors:
RuntimeError: The communication pipe between the main process and its spawned children is broken.
In Windows OS, this usually means that the child process raised an exception while it was being spawned, before it was setup to communicate to the main process.
The exceptions raised by the child process while spawning cannot be caught or handled from the main process, and when running from an IPython or jupyter notebook interactive kernel, the child’s exception and traceback appears to be lost.
A known way to see the child’s error, and try to fix or handle it, is to run the problematic code as a batch script from a system’s Command Prompt. The child’s exception will be printed to the Command Promt’s stderr, and it should be visible above this error and traceback.
Note that if running a jupyter notebook that was invoked from a Command Prompt, the child’s exception should have been printed to the Command Prompt on which the notebook is running.
Could you please check what is incorrect?
import pymc3 as pm
import numpy as np
import scipy.stats as stats
#simulate samples from a bernoulli variable: feat1
p0 = 0.5 #true p of the Bernoulli
obs_feat1 = stats.bernoulli.rvs(0.5, size=1000)
#simulate samples from a normal variable: feat2
mu0 = 0 #true mu of Gaussian
sigma0 = 1 #true sigma of Gaussian
obs_feat2 = stats.norm.rvs(loc=mu0, scale=sigma0, size=1000)
#use Densitydist to get customized 'joint' distribution for the two random varaibles
#logp_ = w1 * ln( LL(feat1) ) + w2 * ln( LL(feat2) )
(w1,w2) = (0.5,0.5)
import theano.tensor as tt
def log_likelihood(p_1, mu_2, sigma_2, f1, f2, w1, w2):
#p_1: parameter for bernoulli variable (1st variable)
#mu_2, sigma_2: parameters for normal variable (2nd varaiable)
#f1: an observation of feat1
#f2: an observation of feat2
#w1, w2: weight for building likelihood
LL1 = tt.switch( f1, tt.log(p_1), tt.log(1 - p_1) )
tau_2 = 1. / (sigma_2**2)
LL2 = (-tau_2 * (f2 - mu_2)**2 + tt.log(tau_2 / np.pi / 2.)) / 2.
logP = sum(w1*LL1 + w2*LL2)
with pm.Model() as model:
p_1 = pm.Uniform('p_1', lower=0, upper=1)
mu_2 = pm.Uniform('mu_2', lower=-1, upper=1)
sigma_2 = pm.Uniform('mu_2', lower=0, upper=2)
like = pm.DensityDist('like', log_likelihood,
w1 = w1,
w2 = w2,)
step = pm.Metropolis()
trace = pm.sample(10000,step=step)