Sampling failure using initial values with Metropolis

[fonnesbeck]

Running certain models using the current master breaks when initial values are passed along with Metropolis as a step method. Works with no initial values, and with NUTS + inits.

from pymc3 import Normal, Binomial, sample, Model, Metropolis
from pymc3.math import exp
import numpy as np

# Samples for each dose level
n = 5 * np.ones(4, dtype=int)
# Log-dose
dose = np.array([-.86, -.3, -.05, .73])
deaths = np.array([0, 1, 3, 5])

def invlogit(x):
    return exp(x) / (1 + exp(x))

with Model() as bioassay_model:

    # Logit-linear model parameters
    alpha = Normal('alpha', 0, 0.01)
    beta = Normal('beta', 0, 0.01)

    # Calculate probabilities of death
    theta = invlogit(alpha + beta * dose)

    # Data likelihood
    deaths = Binomial('deaths', n=n, p=theta, observed=deaths)
    
# This works
with bioassay_model:
    
    bioassay_trace = sample(2000, step=Metropolis())

# This works   
with bioassay_model:
    
    bioassay_trace = sample(200, step=NUTS(), start={'alpha':1})

# This breaks   
with bioassay_model:
    
    bioassay_trace = sample(2000, step=Metropolis(), start={'alpha':1})

yields:

Traceback (most recent call last):
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/theano/compile/function_module.py", line 873, 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 "broken_sample.py", line 39, in <module>
    bioassay_trace = sample(2000, step=Metropolis(), start={'alpha':1})
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/sampling.py", line 175, in sample
    return sample_func(**sample_args)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/sampling.py", line 185, in _sample
    for strace in sampling:
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/tqdm/_tqdm.py", line 830, in __iter__
    for obj in iterable:
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/sampling.py", line 267, in _iter_sample
    point = step.step(point)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/step_methods/compound.py", line 16, in step
    point = method.step(point)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/step_methods/arraystep.py", line 142, in step
    apoint = self.astep(bij.map(point))
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/pymc3/step_methods/metropolis.py", line 124, in astep
    q_new = metrop_select(self.delta_logp(q, q0), q, q0)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/theano/compile/function_module.py", line 886, in __call__
    storage_map=getattr(self.fn, 'storage_map', None))
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/theano/gof/link.py", line 325, in raise_with_op
    reraise(exc_type, exc_value, exc_trace)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/six.py", line 685, in reraise
    raise value.with_traceback(tb)
  File "/Users/fonnescj/anaconda3/envs/pymc/lib/python3.5/site-packages/theano/compile/function_module.py", line 873, 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{Composite{Switch(Identity((GE(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i3) * LE(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i4))), (i5 + Switch(EQ(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i3), i6, (i7 * i8 * scalar_softplus((-Composite{(i0 + (i1 * i2))}(i0, i1, i2))))) + Switch(EQ((i9 - scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2))), i3), i6, (i10 * i11 * scalar_softplus(Composite{(i0 + (i1 * i2))}(i0, i1, i2))))), i6)}}(InplaceDimShuffle{x}.0, inarray1, TensorConstant{[-0.86 -0....05  0.73]}, TensorConstant{(1,) of 0}, TensorConstant{(1,) of 1}, TensorConstant{[ 0.      ...        ]}, TensorConstant{(1,) of -inf}, TensorConstant{(1,) of -1.0}, TensorConstant{[ 0.  1.  3.  5.]}, TensorConstant{(1,) of 1.0}, TensorConstant{(1,) of -1.0}, TensorConstant{[ 5.  4.  2.  0.]})
Toposort index: 4
Inputs types: [TensorType(float64, (True,)), TensorType(float64, (True,)), TensorType(float64, vector), TensorType(int8, (True,)), TensorType(int8, (True,)), TensorType(float64, vector), TensorType(float32, (True,)), TensorType(float64, (True,)), TensorType(float64, vector), TensorType(float64, (True,)), TensorType(float64, (True,)), TensorType(float64, vector)]
Inputs shapes: [(1,), (1,), (4,), (1,), (1,), (4,), (1,), (1,), (4,), (1,), (1,), (4,)]
Inputs strides: [(8,), (8,), (8,), (1,), (1,), (8,), (4,), (8,), (8,), (8,), (8,), (8,)]
Inputs values: [array([1]), array([ 0.0777477]), array([-0.86, -0.3 , -0.05,  0.73]), array([0], dtype=int8), array([1], dtype=int8), array([ 0.        ,  1.60943791,  2.30258509,  0.        ]), array([-inf], dtype=float32), array([-1.]), array([ 0.,  1.,  3.,  5.]), array([ 1.]), array([-1.]), array([ 5.,  4.,  2.,  0.])]
Outputs clients: [[Sum{acc_dtype=float64}(Elemwise{Composite{Switch(Identity((GE(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i3) * LE(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i4))), (i5 + Switch(EQ(scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2)), i3), i6, (i7 * i8 * scalar_softplus((-Composite{(i0 + (i1 * i2))}(i0, i1, i2))))) + Switch(EQ((i9 - scalar_sigmoid(Composite{(i0 + (i1 * i2))}(i0, i1, i2))), i3), i6, (i10 * i11 * scalar_softplus(Composite{(i0 + (i1 * i2))}(i0, i1, i2))))), i6)}}.0)]]

[fonnesbeck]

This is due to an integer being used as the initial value. These should be cast automatically.

Not sure why it breaks for Metropolis but not NUTS. Is NUTS using the starting value?

[twiecki]

Related: #1640
types are a big problem.

[junpenglao]

I just test this on master and seems there is no bug anymore. Closing...