I'd like to model a distribution which is a mixture of a Normal and the constant 0.
I couldn't find a solution because in all the mixture examples I've found the class of distribution is the same for every category.
Here is some code to illustrate what I'm looking for:
with pm.Model() as model:
x_non_zero = pm.Normal(...)
zero_rate = pm.Uniform('zero_rate', lower=0.0, upper=.0, testval=0.5)
fr = pm.Bernoulli('fr', p=zero_rate)
x = pm.???('x', pm.switch(pm.eq(fr, 0), x_non_zero, 0), observed=data['x'])
I'm interested in the rate the data is exactly zero and the parameters of the normal when it is non-zero.
One option will be to try with a Gaussian mixture model, we may think of a Gaussian with sd=0
as a constant value. Another option will be to use a model like the following:
with pm.Model() as model:
mean = pm.Normal('mean', mu=100, sd=10)
sd = pm.HalfNormal('sd', sd=10)
category = pm.Categorical('category', p=[0.5, 0.5], shape=len(x))
mu = pm.switch(pm.eq(category, 0), 0, mean)
eps = pm.switch(pm.eq(category, 0), 0.1, sd)
obs = pm.Normal('obs', mu=mu, sd=eps, observed=x)
step0 = pm.ElemwiseCategorical(vars=[category], values=[0, 1])
step1 = pm.Metropolis(vars=[mean, sd])
trace = pm.sample(10000, step=[step0, step1])
to find out the rate you can compute
burnin = 100
np.mean(trace[burnin]['category'])