TypeError with blocked `Metropolis`

Hello,

I need to Metropolis-update two (very) correlated random variables together. One of the two RVs is discrete, so I can’t use variational methods for it. I think a good strategy is to use blocked=True flag in the step method, but I get the TypeError below.

I can recreate the error in the disaster model example:

disasters_data = array([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
                3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
                2, 2, 3, 4, 2, 1, 3, 2, 2, 1, 1, 1, 1, 3, 0, 0,
                1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
                0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
                3, 3, 1, 1, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
                0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
year = arange(1851, 1962)
with pm.Model() as model:
  switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max())
  early_mean = pm.Exponential('early_mean', lam=1.)
  late_mean = pm.Exponential('late_mean', lam=1.)
  rate = tt.switch(switchpoint >= year, early_mean, late_mean)
  disasters = pm.Poisson('disasters', rate, observed=disasters_data)

  blocked_step = pm.Metropolis([switchpoint, early_mean], blocked=True)
  tr = pm.sample(1000, tune=500, step=[blocked_step])

Here is the output:

Assigned NUTS to late_mean_log__
  0%|          | 0/1500 [00:00<?, ?it/s]
Traceback (most recent call last):
  File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/theano/compile/function_module.py", line 884, in __call__
   self.fn() if output_subset is None else\
TypeError: expected type_num 12 (NPY_FLOAT64) got 7


During handling of the above exception, another exception occurred:

  Traceback (most recent call last):
    File "disaster_model.py", line 47, in <module>
      tr = pm.sample(1000, tune=500, start=start, step=blocked_step)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/sampling.py", line 286, in sample
      return sample_func(**sample_args)[discard:]
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/sampling.py", line 333, in _sample
      for it, strace in enumerate(sampling):
    File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/tqdm/_tqdm.py", line 872, in __iter__
      for obj in iterable:
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/sampling.py", line 431, in _iter_sample
      point, states = step.step(point)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/step_methods/compound.py", line 24, in step
      point, state = method.step(point)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/step_methods/arraystep.py", line 153, in step
      apoint, stats = self.astep(bij.map(point))
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/step_methods/metropolis.py", line 156, in astep
      accept = self.delta_logp(q, q0)
    File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/theano/compile/function_module.py", line 898, in __call__
      storage_map=getattr(self.fn, 'storage_map', None))
    File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/theano/gof/link.py", line 325, in raise_with_op
      reraise(exc_type, exc_value, exc_trace)
    File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/six.py", line 685, in reraise
      raise value.with_traceback(tb)
    File "/home/u1573864/.conda/envs/PyMC3_dev/lib/python3.6/site-packages/theano/compile/function_module.py", line 884, in __call__
      self.fn() if output_subset is None else\
  TypeError: expected type_num 12 (NPY_FLOAT64) got 7
  Apply node that caused the error: Elemwise{Cast{int64}}(Subtensor{int64}.0)
  Toposort index: 20
  Inputs types: [TensorType(float64, scalar)]
  Inputs shapes: [()]
  Inputs strides: [()]
  Inputs values: [array(1906)]
  Outputs clients: [[InplaceDimShuffle{x}(Elemwise{Cast{int64}}.0), Elemwise{Composite{Switch(Identity((GE(i0, i1) * LE(i0, i2))), i3, i4)}}(Elemwise{Cast{int64}}.0, TensorConstant{1851}, TensorConstant{1961}, TensorConstant{-4.709530201312334}, TensorConstant{-inf})]]

  Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
    File "disaster_model.py", line 45, in <module>
      blocked_step = pm.Metropolis([switchpoint, early_mean], blocked=1)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/step_methods/metropolis.py", line 130, in __init__
      self.delta_logp = delta_logp(model.logpt, vars, shared)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/step_methods/metropolis.py", line 479, in delta_logp
      [logp0], inarray0 = pm.join_nonshared_inputs([logp], vars, shared)
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/theanof.py", line 246, in join_nonshared_inputs
      for var, slc, shp, dtyp in ordering.vmap}
    File "/home/u1573864/PyMC3_dev/pymc3/pymc3/theanof.py", line 246, in <dictcomp>
      for var, slc, shp, dtyp in ordering.vmap}

  HINT: Use the Theano flag 'exception_verbosity=high' for a debugprint and storage map footprint of this apply node.

Probably this is caused by having a single block with both integer and float data types…

PS: this topic is specific to the TypeError with blocked metropolis and is different from a previous one

PPS: there is an open issue with the same error here: https://github.com/pymc-devs/pymc3/issues/1681

Hmm, the problem seems to be that you can not do block=True when you have both continuous and discrete variables. So:

with pm.Model() as model:
  switchpoint = pm.Uniform('switchpoint', lower=year.min(), upper=year.max())
  blocked_step = pm.Metropolis([switchpoint, early_mean], blocked=True)

#Works

with pm.Model() as model:
  switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max())
  blocked_step = pm.Metropolis([switchpoint, early_mean], blocked=False)

#Works

with pm.Model() as model:
  switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max())
  blocked_step = pm.Metropolis([switchpoint, early_mean], blocked=True)

#Does not Works

In PyMC2 that was possible. How hard is to code that in PyMC3? Is it OK to open an issue on pymc-devs?

PR welcome :wink: I think it might work if we internally cast the int to float, perform the update, and cast that back to int.

I fear that the problem involves type subleties in theano as in the issue 1681, but it’s worth having a look. I opened a new issue here: https://github.com/pymc-devs/pymc3/issues/2525.