Multinomial dlogpt produces nan - 3.0.rc2 Regression

[asford]

Multinomial.dlogpt results in nan for valid inputs in pymc3 version 3.0.rc2. This is a regression from version 3.0.rc1, and reproduces on both debian and osx. I've attached a minimal reproduction below.

I don't have much experience with debugging pymc3 or theano, but would be happy to provide additional information or submit a reproduction test case if it would be useful.

Reproduction

from __future__ import print_function
import numpy
import theano
import pymc3
from pymc3.theanof import gradient

import os
print("uname", os.uname()[2:])
print("numpy", numpy.__version__)
print("theano", theano.__version__)
print("pymc3", pymc3.__version__)

with pymc3.Model() as dbg_model:
    numpy.random.seed(1663)
    obs = numpy.floor(numpy.random.random(3) * 1000).astype(int)

    infrac = pymc3.Normal("infrac", sd=1, shape=3, testval=.1)
    infrac = infrac / infrac.sum()

    mn = pymc3.Multinomial(
        name = "mn",
        n=obs.sum(),
        p=infrac,
        observed=obs
    )

    dlogpt = dbg_model.fn(gradient(dbg_model.logpt, None))
    logpt = dbg_model.fn(dbg_model.logpt)

    print("logpt", logpt( **dbg_model.test_point ))
    print("dlogpt", dlogpt(**dbg_model.test_point ))

3.0.rc2 Result

uname ('3.16.0-4-amd64', '#1 SMP Debian 3.16.7-ckt25-2 (2016-04-08)', 'x86_64')
numpy 1.11.0
theano 0.8.2
pymc3 3.0.rc2
logpt -25.0047103337
dlogpt [ nan  nan  nan]

3.0.rc1 Expected Result

uname ('3.16.0-4-amd64', '#1 SMP Debian 3.16.7-ckt25-2 (2016-04-08)', 'x86_64')
numpy 1.11.0
theano 0.8.2
pymc3 3.0.rc1
logpt -25.0047103337
dlogpt [ 1133.23333333   -96.76666667 -1036.76666667]

[twiecki]

The only place this could have been introduced is #1390. CC @fonnesbeck

[fonnesbeck]

I will check it out. BTW, using a normal for infrac is a very poor choice, as it will result in proposals that are not greater than zero, making sampling inefficient. I assume you meant to use a half-normal?

[fonnesbeck]

I notice that the testval for infrac is 0.1 for a variable with shape=3. We need to make this raise an exception if it does not. You then clobber infrac by assigning its scaled value in place of the original. I would not expect this model to work, but we need to provide the appropriate error message.

[asford]

This is just a minimal repro case, it's not intended as a meaningful model. Thank you for investigating this.

[asford]

In 3.0.rc1 a scalar testval broadcasts into the appropriately shaped array. That behavior is useful and should be preserved.

from __future__ import print_function
import numpy
import theano
import pymc3
from pymc3.theanof import gradient

import os
print("uname", os.uname()[2:])
print("numpy", numpy.__version__)
print("theano", theano.__version__)
print("pymc3", pymc3.__version__)

with pymc3.Model() as dbg_model:
    numpy.random.seed(1663)
    obs = numpy.floor(numpy.random.random(3) * 1000).astype(int)

    infrac = pymc3.Normal("infrac", sd=1, shape=3, testval=.1)
    infrac = infrac / infrac.sum()

    mn = pymc3.Multinomial(
        name = "mn",
        n=obs.sum(),
        p=infrac,
        observed=obs
    )

    dlogpt = dbg_model.fn(gradient(dbg_model.logpt, None))
    logpt = dbg_model.fn(dbg_model.logpt)

    print("logpt", logpt( **dbg_model.test_point ))
    print("dlogpt", dlogpt(**dbg_model.test_point ))
    print("test_point", dbg_model.test_point)

uname ('15.6.0', 'Darwin Kernel Version 15.6.0: Mon Aug 29 20:21:34 PDT 2016; root:xnu-3248.60.11~1/RELEASE_X86_64', 'x86_64')
numpy 1.11.2
theano 0.8.2
pymc3 3.0.rc1
logpt -25.0047103337
dlogpt [ 1133.23333333   -96.76666667 -1036.76666667]
test_point {'infrac': array([ 0.1,  0.1,  0.1])}

[fonnesbeck]

This result comes from the following bound test that I added to the distribution:

tt.isclose(p.sum(), 1)

Interestingly, this fails only on the gradient calculation. I would think we would want to enforce this condition.

[fonnesbeck]

Please test the pull request. Works here.