Zero inflated binomial #2251
Zero inflated binomial #2251
Conversation
pymc3/distributions/discrete.py
Outdated
| p : float | ||
| Probability of success in each trial (0 < p < 1). | ||
| psi : float | ||
| Expected proportion of Poisson variates (0 < psi < 1) |
| g = generate_samples(stats.binom.rvs, n, p, | ||
| dist_shape=self.shape, | ||
| size=size) | ||
| sampled = g * (np.random.random(np.squeeze(g.shape)) < psi) |
There was a problem hiding this comment.
What is the squeeze about? Doesn't this break broadcasting with g?
There was a problem hiding this comment.
Doesn't appear to. I'm just mirroring what is happening in the other ZI dists.
There was a problem hiding this comment.
It passes the test_distribution_random test.
pymc3/distributions/discrete.py
Outdated
| def logp(self, value): | ||
| return tt.switch(value > 0, | ||
| tt.log(self.psi) + self.bin.logp(value), | ||
| tt.log((1. - self.psi) + self.psi * tt.pow(1 - self.p, self.n))) |
There was a problem hiding this comment.
Does this work numerically for the gradient? Maybe a test that this is reasonable for corner cases for p and psi and large n would be nice.
There was a problem hiding this comment.
Isnt that what the tests do?
There was a problem hiding this comment.
I think this work for the gradient but is numerically unstable. Why not use the log-sum-exp trick as in pm.Mixture? https://github.com/pymc-devs/pymc3/blob/master/pymc3/distributions/mixture.py#L110-L115
There was a problem hiding this comment.
Ah, right. I was just following the implementations of the other ZI distributions, but they should be "robustified", yes.
There was a problem hiding this comment.
Should be good to go now
pymc3/distributions/discrete.py
Outdated
| n = self.n | ||
|
|
||
| logp_val = tt.switch(value > 0, | ||
| logsumexp(tt.log(psi) + self.bin.logp(value)), |
There was a problem hiding this comment.
Shouldn't this still be just tt.log(self.psi) + self.bin.logp(value)?
There was a problem hiding this comment.
Doesn't logsumexp help here?
There was a problem hiding this comment.
The exp is missing :-)
It helps if you have an expression like log(exp(a) + exp(b)). But we are only doing log(a + b) here.
pymc3/distributions/discrete.py
Outdated
|
|
||
| logp_val = tt.switch(value > 0, | ||
| logsumexp(tt.log(psi) + self.bin.logp(value)), | ||
| logsumexp(tt.log((1. - psi) + psi * tt.pow(1 - p, n)))) |
There was a problem hiding this comment.
I think this should be
logaddexp(tt.log1p(-psi), tt.log(psi) + n * tt.log1p(-p))
where
def logaddexp(a, b):
diff = b - a
return tt.switch(diff > 0, b + tt.log1p(tt.exp(-diff)), a + tt.log1p(tt.exp(diff)))There was a problem hiding this comment.
Actually, our logsumexp doesn't have the same signature, but you are right in principle.
There was a problem hiding this comment.
I think we should just add the logaddexp function. This is nicer (and I think faster) if we have only two variables and comes in handy quite often. I'm surprised that I can't find this in theano already....
There was a problem hiding this comment.
In fact, I think logsumexp is used incorrectly in Mixture, if I am reading it correctly.
There was a problem hiding this comment.
logsumexp is getting a weight vector [psi, 1.-psi], maybe that's why?
| """ | ||
|
|
||
| def __init__(self, theta, psi, *args, **kwargs): | ||
| def __init__(self, psi, theta, *args, **kwargs): |
There was a problem hiding this comment.
isn't this a backward incompatible change?
There was a problem hiding this comment.
This harmonizes them with the convention in Mixture, which is where they will eventually end up. So, a break is coming in one place or the other.
There was a problem hiding this comment.
Then we should definitely put that in the release notes. And maybe also print a warning until 3.2?
|
Agree with @aseyboldt, the logp is off, for example: x = np.concatenate([np.random.poisson(4, size=180), np.zeros(20)])
with pm.Model() as model0:
ψ = pm.Beta('ψ', 1., 1.)
θ = pm.Gamma('θ', 1., 1.)
like = ZeroInflatedPoisson('like', psi=ψ, theta=θ, observed=x)
tr0 = pm.sample(3000, init=None, njobs=2)
pm.traceplot(tr0);gives: |
|
OK, getting something more reasonable now:
|
|
LGTM |
|
Tested locally, works perfect (even the case i showed you few days ago in a nb @fonnesbeck) |
Adds zero-inflated binomial distribution, following the other zero-inflated distributions in
discrete.py. Will eventually be superseded by #2246 but that PR is not yet working. This will make the main zero-inflated distributions available for 3.1.