search

pymc-devs / pymc

Bayesian Modeling and Probabilistic Programming in Python
https://docs.pymc.io/
Other
8.64k stars 1.99k forks source link

Broadcastable pattern error using Metropolis in large graph #1983

Closed junpenglao closed 4 years ago

junpenglao commented 7 years ago

I am reworking on my PyMC3 port of Lee and Wagenmakers' Bayesian Cognitive Modeling book, seems there is a bug running Metropolis sampler if the theano graph is large:

import numpy as np
import theano.tensor as tt
import pymc3 as pm

ns = 20
nc = 9
nq = 30
qcc = np.asarray( [[ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  1.,  1., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  0.,  1.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
                   [ 0.,  1.,  0.,  0.,  0., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  0.,  1., -1.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  1.,  0.,  0.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  0.,  0.,  1.,  1.,  0.,  0.,  0.],
                   [ 0.,  0.,  0.,  1.,  1.,  1.,  0.,  1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
                   [ 0.,  0.,  1., -1.,  0., -1.,  0.,  0., -1.],
                   [ 0.,  0.,  0.,  0.,  0.,  1., -1., -1.,  0.],
                   [ 0.,  0.,  0.,  0.,  0.,  1., -1., -1.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.],
                   [ 0.,  0.,  1.,  0., -1., -1.,  0.,  0.,  0.]] )
qccmat = np.tile(qcc[np.newaxis, :, :], (ns, 1, 1))
y = np.asarray([[1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 0, 1, 1, 0, 1, 1, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 0, 0, 0, 1, 0],
               [1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1,
                1, 1, 0, 1, 0, 1, 1, 1],
               [0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0,
                1, 0, 1, 0, 1, 0, 0, 0],
               [1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 0, 1, 1, 1, 0, 0],
               [1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 1, 1, 0, 0],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 0, 1, 1, 1],
               [1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 1, 1,
                1, 1, 0, 0, 1, 1, 1, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 0],
               [0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1,
                1, 1, 1, 0, 1, 0, 0, 0],
               [0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1,
                1, 1, 0, 1, 0, 0, 1, 0],
               [1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1,
                1, 1, 0, 1, 1, 0, 0, 0],
               [1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 0, 1, 0, 0, 0],
               [1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 1,
                1, 1, 0, 0, 1, 0, 0, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
                1, 1, 1, 1, 1, 1, 1, 1],
               [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1,
                1, 1, 0, 1, 1, 1, 1, 1],
               [1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1,
                0, 1, 1, 0, 1, 1, 0, 0],
               [1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 1, 1,
                1, 1, 0, 0, 0, 0, 0, 1],
               [1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1,
                0, 1, 1, 1, 0, 1, 1, 1]], dtype=int)

vtmp = np.random.randn(ns, 1, nc)**2
with pm.Model() as model3:
    gamma = pm.Uniform('gamma', lower=.5, upper=1)
    gammat = tt.stack([1-gamma, .5, gamma])

    v1 = pm.HalfNormal('v1', sd=1, shape=(ns, 1, nc), testval=vtmp)
    s1 = pm.Deterministic('s1', tt.argsort(v1, axis=2))
    smat2 = tt.tile(s1, (1, nq, 1))  # s[1:nc] <- rank(v[1:nc])

    # TTB Model For Each Question
    ttmp = tt.sum(qccmat*tt.power(2, smat2), axis=2)
    tmat = -1*(-ttmp > 0)+(ttmp > 0)+1

    yiq = pm.Bernoulli('yiq', p=gammat[tmat], observed=y)

    trace3 = pm.sample(1e4, step=pm.Metropolis())

This model compiled without problem, and sample without problem in ADVI and NUTS, however, if I set step=pm.Metropolis() the following error returns:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-14-87778767da2a> in <module>()
     14     yiq = pm.Bernoulli('yiq', p=gammat[tmat], observed=y)
     15 
---> 16     trace3 = pm.sample(1e4, step=pm.Metropolis())

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/arraystep.py in __new__(cls, *args, **kwargs)
     61                 # If we don't return the instance we have to manually
     62                 # call __init__
---> 63                 step.__init__([var], *args, **kwargs)
     64                 # Hack for creating the class correctly when unpickling.
     65                 step.__newargs = ([var], ) + args, kwargs

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py in __init__(self, vars, S, proposal_dist, scaling, tune, tune_interval, model, mode, **kwargs)
    122 
    123         shared = pm.make_shared_replacements(vars, model)
--> 124         self.delta_logp = delta_logp(model.logpt, vars, shared)
    125         super(Metropolis, self).__init__(vars, shared)
    126 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py in delta_logp(logp, vars, shared)
    469 
    470 def delta_logp(logp, vars, shared):
--> 471     [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
    472 
    473     tensor_type = inarray0.type

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/theanof.py in join_nonshared_inputs(xs, vars, shared, make_shared)
    236     replace.update(shared)
    237 
--> 238     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    239     return xs_special, inarray
    240 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/pymc3/theanof.py in <listcomp>(.0)
    236     replace.update(shared)
    237 
--> 238     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    239     return xs_special, inarray
    240 

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/scan_module/scan_utils.py in clone(output, replace, strict, share_inputs, copy_inputs)
    256                                         [],
    257                                         strict,
--> 258                                         share_inputs)
    259 
    260     return outs

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
    230     else:
    231         if isinstance(outputs, Variable):
--> 232             cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
    233             cloned_outputs = cloned_v
    234             # computed_list.append(cloned_v)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(
---> 96                     [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
     97                 for old_o, new_o in zip(owner.outputs, clone_d[owner].outputs):
     98                     clone_d.setdefault(old_o, new_o)

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/gof/graph.py in clone_with_new_inputs(self, inputs, strict)
    240                     remake_node = True
    241         if remake_node:
--> 242             new_node = self.op.make_node(*new_inputs)
    243             new_node.tag = copy(self.tag).__update__(new_node.tag)
    244         else:

/Library/Frameworks/Python.framework/Versions/3.5/lib/python3.5/site-packages/theano/tensor/elemwise.py in make_node(self, _input)
    201                         "The broadcastable pattern of the "
    202                         "input is incorrect for this op. Expected %s, got %s."
--> 203                         % (self.input_broadcastable, ib)))
    204                 # else, expected == b or expected is False and b is True
    205                 # Both case are good.

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (True, False, True, False, True, False), got (True, False, True, False, False, False).

I am on the master branch of pymc3. The background of the model could be found in Chapter 18 of the book.

junpenglao commented 7 years ago

Seems to similar to #1304

junpenglao commented 7 years ago

Found a solution to this problem.

instead of doing:

v1 = pm.HalfNormal('v1', sd=1, shape=(ns, 1, nc))
s1 = pm.Deterministic('s1', tt.argsort(v1, axis=2))

create a variable and then reshape it resolved the problem

v1 = pm.HalfNormal('v1', sd=1, shape=ns*nc)
s1 = tt.argsort(v1.reshape((ns, 1, nc)), axis=2)
twiecki commented 7 years ago

Glad it's working. I don't see why the above code shouldn't work. It could point to a subtle problem in join_nonshared_inputs for higer dimensional inputs, so I'll reopen this.

lerkoah commented 7 years ago

Hi, I have the same problem for my model. 

with rbf_model:
    alpha_model = pm.Normal('alpha', mu=init_alpha, sd=200, shape = numberOfBasis)
    Cx_model = pm.Normal('Cx', mu=init_Cx, sd=20, shape = numberOfBasis)
    Cy_model = pm.Normal('Cy', mu=init_Cy, sd=20, shape = numberOfBasis)
    l_model = pm.Lognormal('l', mu = l, sd=20, shape = 1)

    sigma_model = pm.HalfNormal('sigma', tau = np.ones(2), shape=(2,2), testval=init_sigma)

    PHI_Re, PHI_Im = Vobs_function(U, l_model, alpha_model, Cx_model, Cy_model)
    V_model = tt.stack([PHI_Re, PHI_Im], axis = 1)

    V_obs = pm.MvNormal('V_obs', mu=V_model, cov=sigma_model, observed=V)

    n_samples = 5000
    step = pm.Metropolis()
    trace = pm.sample(n_samples, step)

In my case, my model V_obs is a two dimensional vector, i.e., each measurement has two components. The function Vobs_function() returns this two components and then I concatenate (stack).

I don't know where is the shape-error. Can you help me?

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-63-3118a810fbc9> in <module>()
     47 
     48     n_samples = 5000
---> 49     step = pm.Metropolis()
     50     trace = pm.sample(n_samples, step)
     51 #     estimation = pm.find_MAP()

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py in __new__(cls, *args, **kwargs)
     58                 # If we don't return the instance we have to manually
     59                 # call __init__
---> 60                 step.__init__([var], *args, **kwargs)
     61                 # Hack for creating the class correctly when unpickling.
     62                 step.__newargs = ([var], ) + args, kwargs

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in __init__(self, vars, S, proposal_dist, scaling, tune, tune_interval, model, **kwargs)
     96 
     97         shared = pm.make_shared_replacements(vars, model)
---> 98         self.delta_logp = delta_logp(model.logpt, vars, shared)
     99         super(Metropolis, self).__init__(vars, shared)
    100 

/root/anaconda3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py in delta_logp(logp, vars, shared)
    424 
    425 def delta_logp(logp, vars, shared):
--> 426     [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
    427 
    428     tensor_type = inarray0.type

/root/anaconda3/lib/python3.6/site-packages/pymc3/theanof.py in join_nonshared_inputs(xs, vars, shared, make_shared)
    198     replace.update(shared)
    199 
--> 200     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    201     return xs_special, inarray
    202 

/root/anaconda3/lib/python3.6/site-packages/pymc3/theanof.py in <listcomp>(.0)
    198     replace.update(shared)
    199 
--> 200     xs_special = [theano.clone(x, replace, strict=False) for x in xs]
    201     return xs_special, inarray
    202 

/root/anaconda3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py in clone(output, replace, strict, share_inputs, copy_inputs)
    256                                         [],
    257                                         strict,
--> 258                                         share_inputs)
    259 
    260     return outs

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in rebuild_collect_shared(outputs, inputs, replace, updates, rebuild_strict, copy_inputs_over, no_default_updates)
    230     else:
    231         if isinstance(outputs, Variable):
--> 232             cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over)
    233             cloned_outputs = cloned_v
    234             # computed_list.append(cloned_v)

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     91             if owner not in clone_d:
     92                 for i in owner.inputs:
---> 93                     clone_v_get_shared_updates(i, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(

/root/anaconda3/lib/python3.6/site-packages/theano/compile/pfunc.py in clone_v_get_shared_updates(v, copy_inputs_over)
     94 
     95                 clone_d[owner] = owner.clone_with_new_inputs(
---> 96                     [clone_d[i] for i in owner.inputs], strict=rebuild_strict)
     97                 for old_o, new_o in zip(owner.outputs, clone_d[owner].outputs):
     98                     clone_d.setdefault(old_o, new_o)

/root/anaconda3/lib/python3.6/site-packages/theano/gof/graph.py in clone_with_new_inputs(self, inputs, strict)
    240                     remake_node = True
    241         if remake_node:
--> 242             new_node = self.op.make_node(*new_inputs)
    243             new_node.tag = copy(self.tag).__update__(new_node.tag)
    244         else:

/root/anaconda3/lib/python3.6/site-packages/theano/tensor/elemwise.py in make_node(self, _input)
    201                         "The broadcastable pattern of the "
    202                         "input is incorrect for this op. Expected %s, got %s."
--> 203                         % (self.input_broadcastable, ib)))
    204                 # else, expected == b or expected is False and b is True
    205                 # Both case are good.

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (True,), got (False,).
junpenglao commented 7 years ago

hmm, it is a bit difficult to say without the data and the parameters. But I dont think you should use sigma_model = pm.HalfNormal('sigma', tau = np.ones(2), shape=(2,2), testval=init_sigma) as your prior for cov. Try pm.LKJCholeskyCov it should work better

lerkoah commented 7 years ago

Using LKJCholeskyCov the problem doesn't disappear.

I think that is a problem from the source specificly for the function join_nonshared_inputs. I'm going to make a fork for trying to solve this issue. For me, It is a persona challenge ahha.

giovannidoni commented 7 years ago

I am the same problem here when trying to do some hierarchical stuff.

ids = D.index.values.astype(int)

X = D.values[:,1:-2]

with pm.Model() as maxdiff:

    mu = pm.Normal('hyper_alpha', mu=0, sd=5, shape=1)
    sigma = pm.HalfCauchy('hyper_sigma', beta=5)  

    alpha = pm.Normal('alpha', mu=mu, sd=sigma, shape=(dd.ids, dd.attr)) 
    coeff = pm.Normal('beta', mu=mu, sd=sigma, shape=(dd.ids, dd.attr)) 

    theta = coeff[ids] * X + alpha[ids]

    p = pm.Deterministic('p', T.nnet.softmax(theta))

    y_b = pm.Categorical('yl', p=p, observed=D['best'].values)     

    step = pm.Metropolis(vars=[alpha, coeff])
    trace = pm.sample(50000, step)

Mocked input have right dimensions:

X = D.values[:,1:-2]
print X.shape

coeff = np.random.normal(size=(dd.ids, dd.attr))
coeff = coeff[ids]
print coeff.shape

print (coeff * X).shape

Output: 

(2400, 12)
(2400, 12)
(2400, 12)

The error traceback is the same as previous comments or at least the error reads analogously. The non hierarchical model works fine. Is it a problem generated by the indexing (which are actually used also in previous cases)?

mfansler commented 6 years ago

Same issue (crossposting on SO). It's not clear to me how the workaround for other models in this thread translates to my model implementation. If anyone can help with that, I'd greatly appreciate it!

Model

G = 10    # genes
K = 5     # cell types
C_k = 20  # cells per type

data = sample_data(G, K, C_k)

with pm.Model() as capture_efficiency:

    # Genes expression levels per cell type
    mu_gk = pm.Uniform('mu_gk', 1, 1000, shape=(G, K, 1))

    # Cell capture efficiencies
    p_kc = pm.Beta('p_kc', shape=(1, K, C_k), alpha=5, beta=95)

    # Captured transcripts
    N_gkc = pm.Binomial('N_gkc', shape=(G, K, C_k),
                        n=tt.tensordot(mu_gk, np.ones((C_k, 1)), [[2], [1]]),
                        p=tt.tensordot(np.ones((G, 1)), p_kc, [[1], [0]]),
                        observed=data)

    trace = pm.sample(5000, tune=10000, step=pm.Metropolis())

Sample Data

def sample_data(G=1, K=2, C_k=100):
    mu_gk = np.random.randint(1, 1000, size=(G, K))
    p_kc = np.random.beta(5, 95, (K, C_k))
    N_gkc = np.random.binomial(mu_gk[:, :, np.newaxis], p_kc[np.newaxis, :, :])

    return N_gkc

Error and Stack Trace

Traceback (most recent call last): File "/Applications/PyCharm.app/Contents/helpers/pydev/pydev_run_in_console.py", line 52, in run_file pydev_imports.execfile(file, globals, locals) # execute the script File "/Applications/PyCharm.app/Contents/helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile exec(compile(contents+"\n", file, 'exec'), glob, loc) File "/Users/mfansler/Projects/pymc3/intro/capture-efficiency-celltypes.py", line 46, in trace = pm.sample(5000, tune=10000, step=pm.Metropolis(), File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 65, in new step.__init__([var], *args, **kwargs) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 136, in init self.delta_logp = delta_logp(model.logpt, vars, shared) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/step_methods/metropolis.py", line 624, in delta_logp [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in join_nonshared_inputs xs_special = [theano.clone(x, replace, strict=False) for x in xs] File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/pymc3/theanof.py", line 245, in xs_special = [theano.clone(x, replace, strict=False) for x in xs] File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/scan_module/scan_utils.py", line 247, in clone share_inputs) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 232, in rebuild_collect_shared cloned_v = clone_v_get_shared_updates(outputs, copy_inputs_over) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates clone_v_get_shared_updates(i, copy_inputs_over) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates clone_v_get_shared_updates(i, copy_inputs_over) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 93, in clone_v_get_shared_updates clone_v_get_shared_updates(i, copy_inputs_over) [Previous line repeated 9 more times] File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/compile/pfunc.py", line 96, in clone_v_get_shared_updates [clone_d[i] for i in owner.inputs], strict=rebuild_strict) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/gof/graph.py", line 246, in clone_with_new_inputs new_node = self.op.make_node(*new_inputs) File "/Users/mfansler/anaconda/envs/pymc3/lib/python3.6/site-packages/theano/tensor/elemwise.py", line 230, in make_node % (self.input_broadcastable, ib)))

TypeError: The broadcastable pattern of the input is incorrect for this op. Expected (False, False, True), got (False, False, False).

natalieklein commented 5 years ago

I know this thread is older, but I just had the same issue and wanted to share what worked for me. Oddly, what fixed it was not flattening anything (tried that, didn't fix it). It turns out I had one variable that had shape = 1 explicitly specified (because it used to have a different shape, otherwise I normally wouldn't explicitly specify it). Taking that shape = 1 away seems to fix Metropolis sampling.

porterjenkins commented 5 years ago

@natalieklein Thanks for posting! That seemed to work for me as well. Changing my shapes from (N x 1) --> (N) seemed to fix the issue for me.

twiecki commented 4 years ago

Feel free to reopen if still an issue.