Erroneous dimension mis-match error

[fonnesbeck]

As reported by another user, the following simple model fails on OS X and Python 3. It appears to be complaining about the size of the observed argument for the NegativeBinomial, which should not be happening. Changing to another distribution (e.g. Poisson) does not change this. Adding a shape argument to the likelihood also does not fix the problem, nor does removing find_map. Passing a scalar as data resolves the issue.

import pymc3 as pm
import numpy as np

data = [6.0,12.0,12.0,46.0,5.0,11.0,11.0,39.0,4.0,10.0,25.0,11.0,8.0,5.0,10.0,2.0,30.0,21.0]

with pm.Model() as model:
    alpha = pm.Uniform('alpha', lower=0.1, upper=100)
    mu = pm.Uniform('mu', lower=0.1, upper=100)

    y_pred = pm.NegativeBinomial('y_pred', mu=mu, alpha=alpha)

    y_est = pm.NegativeBinomial('y_est', 
                                mu=mu, 
                                alpha=alpha, 
                                observed=data)

    start = pm.find_MAP()

    step = pm.NUTS()
    trace = pm.sample(2000, step, start=start, progressbar=True)

Here's the error.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
/Users/fonnescj/anaconda3/lib/python3.4/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    839         try:
--> 840             outputs = self.fn()
    841         except Exception:

ValueError: Input dimension mis-match. (input[0].shape[0] = 1, input[4].shape[0] = 18)

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-4-899599b91d18> in <module>()
     15                                 observed=data, shape=len(data))
     16 
---> 17     start = pm.find_MAP()
     18 
     19     step = pm.NUTS()

/Users/fonnescj/Github/pymc3/pymc3/tuning/starting.py in find_MAP(start, vars, fmin, return_raw, disp, model, *args, **kwargs)
     79     if 'fprime' in getargspec(fmin).args:
     80         r = fmin(logp_o, bij.map(
---> 81             start), fprime=grad_logp_o, disp=disp, *args, **kwargs)
     82     else:
     83         r = fmin(logp_o, bij.map(start), disp=disp, *args, **kwargs)

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/scipy/optimize/optimize.py in fmin_bfgs(f, x0, fprime, args, gtol, norm, epsilon, maxiter, full_output, disp, retall, callback)
    787             'return_all': retall}
    788 
--> 789     res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)
    790 
    791     if full_output:

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/scipy/optimize/optimize.py in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
    841     else:
    842         grad_calls, myfprime = wrap_function(fprime, args)
--> 843     gfk = myfprime(x0)
    844     k = 0
    845     N = len(x0)

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/scipy/optimize/optimize.py in function_wrapper(*wrapper_args)
    283     def function_wrapper(*wrapper_args):
    284         ncalls[0] += 1
--> 285         return function(*(wrapper_args + args))
    286 
    287     return ncalls, function_wrapper

/Users/fonnescj/Github/pymc3/pymc3/tuning/starting.py in grad_logp_o(point)
     74 
     75     def grad_logp_o(point):
---> 76         return nan_to_num(-dlogp(point))
     77 
     78     # Check to see if minimization function actually uses the gradient

/Users/fonnescj/Github/pymc3/pymc3/blocking.py in __call__(self, x)
    117 
    118     def __call__(self, x):
--> 119         return self.fa(self.fb(x))

/Users/fonnescj/Github/pymc3/pymc3/model.py in __call__(self, state)
    397 
    398     def __call__(self, state):
--> 399         return self.f(**state)
    400 
    401 class LoosePointFunc(object):

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    849                         self.fn.nodes[self.fn.position_of_error],
    850                         self.fn.thunks[self.fn.position_of_error],
--> 851                         storage_map=self.fn.storage_map)
    852                 else:
    853                     # For the c linker We don't have access from

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/theano/gof/link.py in raise_with_op(node, thunk, exc_info, storage_map)
    312         # extra long error message in that case.
    313         pass
--> 314     reraise(exc_type, exc_value, exc_trace)
    315 
    316 

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/six.py in reraise(tp, value, tb)
    656             value = tp()
    657         if value.__traceback__ is not tb:
--> 658             raise value.with_traceback(tb)
    659         raise value
    660 

/Users/fonnescj/anaconda3/lib/python3.4/site-packages/theano/compile/function_module.py in __call__(self, *args, **kwargs)
    838         t0_fn = time.time()
    839         try:
--> 840             outputs = self.fn()
    841         except Exception:
    842             if hasattr(self.fn, 'position_of_error'):

ValueError: Input dimension mis-match. (input[0].shape[0] = 1, input[4].shape[0] = 18)
Apply node that caused the error: Elemwise{Composite{Switch(i0, i1, Switch(i2, Switch(i3, i1, i4), i1))}}(TensorConstant{(1,) of 0}, TensorConstant{(1,) of 0}, Elemwise{mul,no_inplace}.0, InplaceDimShuffle{x}.0, TensorConstant{[  6.  12... 30.  21.]})
Toposort index: 36
Inputs types: [TensorType(int8, vector), TensorType(int8, (True,)), TensorType(int8, (True,)), TensorType(int8, (True,)), TensorType(float64, vector)]
Inputs shapes: [(1,), (1,), (1,), (1,), (18,)]
Inputs strides: [(1,), (1,), (1,), (1,), (8,)]
Inputs values: [array([0], dtype=int8), array([0], dtype=int8), array([1], dtype=int8), array([0], dtype=int8), 'not shown']
Outputs clients: [[Sum{acc_dtype=float64}(Elemwise{Composite{Switch(i0, i1, Switch(i2, Switch(i3, i1, i4), i1))}}.0)]]

HINT: Re-running with most Theano optimization disabled could give you a back-trace of when this node was created. This can be done with by setting the Theano flag 'optimizer=fast_compile'. If that does not work, Theano optimizations can be disabled with 'optimizer=None'.
HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

[markdregan]

I'm not sure if this is useful whilst debugging/investigating - but if the observed dataset contains the value zero, the model executes without error.

data = [0.0, 6.0,12.0,12.0,46.0,5.0,11.0,11.0,39.0,4.0,10.0,25.0,11.0,8.0,5.0,10.0,2.0,30.0,21.0]

[jsalvatier]

wtf. This is all real weird. Especially the zero thing.

On Thu, Oct 1, 2015 at 8:06 AM Mark Regan notifications@github.com wrote:

I'm not sure if this is useful whilst debugging/investigating - but if the
observed dataset contains the value zero, the mode executes without error.

data = [0.0, 6.0,12.0,12.0,46.0,5.0,11.0,11.0,39.0,4.0,10.0,25.0,11.0,8.0,5.0,10.0,2.0,30.0,21.0]

—
Reply to this email directly or view it on GitHub
#842 (comment).

[keryil]

I have just run into a similar issue on OS X, Python 2.7.

model = Model()
obs = np.array([16, 18, 22, 25, 27])
with model:
    theta = pymc3.Uniform('theta', lower=0, upper=1)
    n = pymc3.DiscreteUniform('n', lower=1, upper=500)
    k = pymc3.Binomial("k", n=n, p=theta, shape=1, 
                       observed=obs)
    start = find_MAP()
    step = pymc3.NUTS()
    trace = sample(1000, step, start=start, 
                  progressbar = True)

Produces a similar error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-51-0647933798bd> in <module>()
      6     k = pymc3.Binomial("k", n=n, p=theta, shape=1, 
      7                        observed=obs)
----> 8     start = find_MAP()
      9     step = pymc3.NUTS()
     10     trace = sample(1000, step, start=start, 

/usr/local/lib/python2.7/site-packages/pymc3/tuning/starting.pyc in find_MAP(start, vars, fmin, return_raw, disp, model, *args, **kwargs)
     79     if 'fprime' in getargspec(fmin).args:
     80         r = fmin(logp_o, bij.map(
---> 81             start), fprime=grad_logp_o, disp=disp, *args, **kwargs)
     82     else:
     83         r = fmin(logp_o, bij.map(start), disp=disp, *args, **kwargs)

/usr/local/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in fmin_bfgs(f, x0, fprime, args, gtol, norm, epsilon, maxiter, full_output, disp, retall, callback)
    775             'return_all': retall}
    776 
--> 777     res = _minimize_bfgs(f, x0, args, fprime, callback=callback, **opts)
    778 
    779     if full_output:

/usr/local/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
    830     else:
    831         grad_calls, myfprime = wrap_function(fprime, args)
--> 832     gfk = myfprime(x0)
    833     k = 0
    834     N = len(x0)

/usr/local/lib/python2.7/site-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    279     def function_wrapper(*wrapper_args):
    280         ncalls[0] += 1
--> 281         return function(*(wrapper_args + args))
    282 
    283     return ncalls, function_wrapper

/usr/local/lib/python2.7/site-packages/pymc3/tuning/starting.pyc in grad_logp_o(point)
     74 
     75     def grad_logp_o(point):
---> 76         return nan_to_num(-dlogp(point))
     77 
     78     # Check to see if minimization function actually uses the gradient

/usr/local/lib/python2.7/site-packages/pymc3/blocking.pyc in __call__(self, x)
    117 
    118     def __call__(self, x):
--> 119         return self.fa(self.fb(x))

/usr/local/lib/python2.7/site-packages/pymc3/model.pyc in __call__(self, state)
    397 
    398     def __call__(self, state):
--> 399         return self.f(**state)
    400 
    401 class LoosePointFunc(object):

/usr/local/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    862                     node=self.fn.nodes[self.fn.position_of_error],
    863                     thunk=thunk,
--> 864                     storage_map=getattr(self.fn, 'storage_map', None))
    865             else:
    866                 # old-style linkers raise their own exceptions

/usr/local/lib/python2.7/site-packages/theano/gof/link.pyc in raise_with_op(node, thunk, exc_info, storage_map)
    312         # extra long error message in that case.
    313         pass
--> 314     reraise(exc_type, exc_value, exc_trace)
    315 
    316 

/usr/local/lib/python2.7/site-packages/theano/compile/function_module.pyc in __call__(self, *args, **kwargs)
    850         t0_fn = time.time()
    851         try:
--> 852             outputs = self.fn()
    853         except Exception:
    854             if hasattr(self.fn, 'position_of_error'):

ValueError: Input dimension mis-match. (input[0].shape[0] = 1, input[2].shape[0] = 5)
Apply node that caused the error: Elemwise{Composite{Switch(i0, i1, Switch(i2, i3, i1))}}(TensorConstant{(1,) of 0}, TensorConstant{(1,) of 0}, Elemwise{Composite{(LE(i0, i1) * i2 * i3)}}.0, TensorConstant{[ 16.  18... 25.  27.]})
Toposort index: 16
Inputs types: [TensorType(int8, vector), TensorType(int8, (True,)), TensorType(int8, vector), TensorType(float64, vector)]
Inputs shapes: [(1,), (1,), (5,), (5,)]
Inputs strides: [(1,), (1,), (1,), (8,)]
Inputs values: [array([0], dtype=int8), array([0], dtype=int8), array([1, 1, 1, 1, 1], dtype=int8), array([ 16.,  18.,  22.,  25.,  27.])]
Outputs clients: [[Sum{acc_dtype=float64}(Elemwise{Composite{Switch(i0, i1, Switch(i2, i3, i1))}}.0)]]

HINT: Re-running with most Theano optimization disabled could give you a back-trace of when this node was created. This can be done with by setting the Theano flag 'optimizer=fast_compile'. If that does not work, Theano optimizations can be disabled with 'optimizer=None'.
HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

Moreover, the error DOES mysteriously go away if I change one of the observed values to 0! I suspect this is the same issue.

[markdregan]

The above issue now returns a different error. This might be useful for debugging.

ValueError: Cannot construct a ufunc with more than 32 operands (requested number were: inputs = 32 and outputs = 1)

This is similar to another bug referenced recently: #850

Versions:
Python 3.5.0 (v3.5.0:374f501f4567, Sep 12 2015, 11:00:19)
theano.version == '0.7.0.dev-RELEASE'
pymc3.version == '3.0'

[markdregan]

Update: there are two issues being referred to here.

1. Cannot construct a ufunc with more than 32 variables #850 has resolved the "Cannot construct a ufunc with more than 32 operands" issue
2. After updating pymc3 and dependancies, the original issue returns: "Input dimension mis-match..."

An ipython notebook with this issue can be found here: https://github.com/markdregan/Bayesian-Modelling-in-Python/blob/master/Section%203.%20Hierarchical%20modelling.ipynb

[fonnesbeck]

How about if you cast your observed data to a Numpy array? I've seen lists cause trouble before, but not for awhile.

[markdregan]

Casting to numpy array doesn't resolve the issue. The error trace is the same.

[twiecki]

This is odd indeed, I suspect a bug in NegativeBinomial.

Here is the minimal code to replicate stand-alone:

import pymc3 as pm
import numpy as np

with pm.Model() as model:
    alpha = pm.Uniform('alpha', lower=0, upper=100)
    mu = pm.Uniform('mu', lower=0, upper=100)

    data = [34.0, 27.0, 3.0, 17.0, 17.0, 29.0, 53.0]
    y_est = pm.NegativeBinomial('y_est', mu=mu, alpha=alpha, observed=np.asarray(data, dtype='int16'))

    start = pm.find_MAP()

Btw. @markdregan your project looks very interesting.

[fonnesbeck]

The NegativeBinomial likelihood needs to change in any case. It currently calculates two likelihoods every time it is called and returns one of them, which is not performant. I will have a look.

[fonnesbeck]

Actually, looks like only one expression gets evaluated by Theano, so disregard the above.

[fonnesbeck]

I can confirm that this also occurs with the Poisson.

[fonnesbeck]

The problem is with the switch statement in logpow. I think I have dealt with it. Please test.

[twiecki]

I tested with the original NB of @markdregan and everything worked well.

[markdregan]

Yes - can confirm issue is resolved. Thought I responded here. Thank you all for digging in to the issue and resolving.