pymc3: Disaster example with deterministic switchpoint function

Question

I'm trying to reproduce coal mining example with deterministic function for switchpoint instead of using theano's switch function. Code:

%matplotlib inline
import matplotlib.pyplot as plt
import pymc3
import numpy as np
import theano.tensor as t
import theano

data = np.hstack((np.random.poisson(15,1000),np.random.poisson(2,100)))
plt.plot(data)

@theano.compile.ops.as_op(itypes=[t.lscalar, t.dscalar,t.dscalar],otypes=[t.dvector])
def rate1(sw,mu1,mu2):
    n = len(data)
    out = np.empty(n)
    out[:sw] = mu1
    out[sw:] = mu2
    return out



with pymc3.Model() as dis:
    switchpoint = pymc3.DiscreteUniform('switchpoint',lower=0, upper=len(data)-1)
    mu1 = pymc3.Exponential('mu1', lam=1.)
    mu2 = pymc3.Exponential('mu2',lam=1.)
    disasters=pymc3.Poisson('disasters', mu=rate1, observed = data)

But this code rise an error:

--------------------------------------------------------------------------- KeyError Traceback (most recent call last) c:\program files\git\theano\theano\tensor\type.py in dtype_specs(self) 266 'complex64': (complex, 'theano_complex64', 'NPY_COMPLEX64') --> 267 }[self.dtype] 268 except KeyError:

KeyError: 'object'

During handling of the above exception, another exception occurred:

TypeError Traceback (most recent call last) c:\program files\git\theano\theano\tensor\basic.py in constant_or_value(x, rtype, name, ndim, dtype) 407 rval = rtype( --> 408 TensorType(dtype=x_.dtype, broadcastable=bcastable), 409 x_.copy(),

c:\program files\git\theano\theano\tensor\type.py in init(self, dtype, broadcastable, name, sparse_grad) 49 self.broadcastable = tuple(bool(b) for b in broadcastable) ---> 50 self.dtype_specs() # error checking is done there 51 self.name = name

c:\program files\git\theano\theano\tensor\type.py in dtype_specs(self) 269 raise TypeError("Unsupported dtype for %s: %s" --> 270 % (self.class.name, self.dtype)) 271

TypeError: Unsupported dtype for TensorType: object

During handling of the above exception, another exception occurred:

TypeError Traceback (most recent call last) c:\program files\git\theano\theano\tensor\basic.py in as_tensor_variable(x, name, ndim) 201 try: --> 202 return constant(x, name=name, ndim=ndim) 203 except TypeError:

c:\program files\git\theano\theano\tensor\basic.py in constant(x, name, ndim, dtype) 421 ret = constant_or_value(x, rtype=TensorConstant, name=name, ndim=ndim, --> 422 dtype=dtype) 423

c:\program files\git\theano\theano\tensor\basic.py in constant_or_value(x, rtype, name, ndim, dtype) 416 except Exception: --> 417 raise TypeError("Could not convert %s to TensorType" % x, type(x)) 418

TypeError: ('Could not convert FromFunctionOp{rate1} to TensorType', )

During handling of the above exception, another exception occurred:

AsTensorError Traceback (most recent call last) in () 14 mu2 = pymc3.Exponential('mu2',lam=1.) 15 #rate1 = pymc3.switch(switchpoint >= np.arange(len(data)), mu1,mu2) ---> 16 disasters=pymc3.Poisson('disasters', mu=rate1, observed = data)

C:\Users\User\Anaconda3\lib\site-packages\pymc3\distributions\distribution.py in new(cls, name, *args, **kwargs) 19 if isinstance(name, str): 20 data = kwargs.pop('observed', None) ---> 21 dist = cls.dist(*args, **kwargs) 22 return model.Var(name, dist, data) 23 elif name is None:

C:\Users\User\Anaconda3\lib\site-packages\pymc3\distributions\distribution.py in dist(cls, *args, **kwargs) 32 def dist(cls, *args, **kwargs): 33 dist = object.new(cls) ---> 34 dist.init(*args, **kwargs) 35 return dist 36

C:\Users\User\Anaconda3\lib\site-packages\pymc3\distributions\discrete.py in init(self, mu, *args, **kwargs) 185 super(Poisson, self).init(*args, **kwargs) 186 self.mu = mu --> 187 self.mode = floor(mu).astype('int32') 188 189 def random(self, point=None, size=None, repeat=None):

c:\program files\git\theano\theano\gof\op.py in call(self, *inputs, **kwargs) 598 """ 599 return_list = kwargs.pop('return_list', False) --> 600 node = self.make_node(*inputs, **kwargs) 601 602 if config.compute_test_value != 'off':

c:\program files\git\theano\theano\tensor\elemwise.py in make_node(self, *inputs) 540 using DimShuffle. 541 """ --> 542 inputs = list(map(as_tensor_variable, inputs)) 543 shadow = self.scalar_op.make_node( 544 *[get_scalar_type(dtype=i.type.dtype).make_variable()

c:\program files\git\theano\theano\tensor\basic.py in as_tensor_variable(x, name, ndim) 206 except Exception: 207 str_x = repr(x) --> 208 raise AsTensorError("Cannot convert %s to TensorType" % str_x, type(x)) 209 210 # this has a different name, because _as_tensor_variable is the

AsTensorError: ('Cannot convert FromFunctionOp{rate1} to TensorType', )

How i handle this?

The second thing - when i'm using the pymc3.switch function like this:

with pymc3.Model() as dis:
    switchpoint = pymc3.DiscreteUniform('switchpoint',lower=0, upper=len(data)-1)
    mu1 = pymc3.Exponential('mu1', lam=1.)
    mu2 = pymc3.Exponential('mu2',lam=1.)

    rate1 = pymc3.switch(switchpoint >= np.arange(len(data)), mu1,mu2)

    disasters=pymc3.Poisson('disasters', mu=rate1, observed = data) 

And next try to sample:

with dis:
    step1 = pymc3.NUTS([mu1, mu2])
    step2 = pymc3.Metropolis([switchpoint])
    trace = pymc3.sample(10000, step = [step1,step2])

I get an error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
c:\program files\git\theano\theano\compile\function_module.py in __call__(self, *args, **kwargs)
    858         try:
--> 859             outputs = self.fn()
    860         except Exception:

TypeError: expected type_num 9 (NPY_INT64) got 7

During handling of the above exception, another exception occurred:

TypeError                                 Traceback (most recent call last)
<ipython-input-4-3247d908f897> in <module>()
      2     step1 = pymc3.NUTS([mu1, mu2])
      3     step2 = pymc3.Metropolis([switchpoint])
----> 4     trace = pymc3.sample(10000, step = [step1,step2])

C:\Users\User\Anaconda3\lib\site-packages\pymc3\sampling.py in sample(draws, step, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    153         sample_args = [draws, step, start, trace, chain,
    154                        tune, progressbar, model, random_seed]
--> 155     return sample_func(*sample_args)
    156 
    157 

C:\Users\User\Anaconda3\lib\site-packages\pymc3\sampling.py in _sample(draws, step, start, trace, chain, tune, progressbar, model, random_seed)
    162     progress = progress_bar(draws)
    163     try:
--> 164         for i, strace in enumerate(sampling):
    165             if progressbar:
    166                 progress.update(i)

C:\Users\User\Anaconda3\lib\site-packages\pymc3\sampling.py in _iter_sample(draws, step, start, trace, chain, tune, model, random_seed)
    244         if i == tune:
    245             step = stop_tuning(step)
--> 246         point = step.step(point)
    247         strace.record(point)
    248         yield strace

C:\Users\User\Anaconda3\lib\site-packages\pymc3\step_methods\compound.py in step(self, point)
     11     def step(self, point):
     12         for method in self.methods:
---> 13             point = method.step(point)
     14         return point

C:\Users\User\Anaconda3\lib\site-packages\pymc3\step_methods\arraystep.py in step(self, point)
    116         bij = DictToArrayBijection(self.ordering, point)
    117 
--> 118         apoint = self.astep(bij.map(point))
    119         return bij.rmap(apoint)
    120 

C:\Users\User\Anaconda3\lib\site-packages\pymc3\step_methods\metropolis.py in astep(self, q0)
    123 
    124 
--> 125         q_new = metrop_select(self.delta_logp(q,q0), q, q0)
    126 
    127         if q_new is q:

c:\program files\git\theano\theano\compile\function_module.py in __call__(self, *args, **kwargs)
    869                     node=self.fn.nodes[self.fn.position_of_error],
    870                     thunk=thunk,
--> 871                     storage_map=getattr(self.fn, 'storage_map', None))
    872             else:
    873                 # old-style linkers raise their own exceptions

c:\program files\git\theano\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 

C:\Users\User\Anaconda3\lib\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 

c:\program files\git\theano\theano\compile\function_module.py in __call__(self, *args, **kwargs)
    857         t0_fn = time.time()
    858         try:
--> 859             outputs = self.fn()
    860         except Exception:
    861             if hasattr(self.fn, 'position_of_error'):

TypeError: expected type_num 9 (NPY_INT64) got 7
Apply node that caused the error: Elemwise{Composite{Switch(GE(i0, i1), i2, i3)}}(InplaceDimShuffle{x}.0, TensorConstant{[   0    1..1098 1099]}, InplaceDimShuffle{x}.0, InplaceDimShuffle{x}.0)
Toposort index: 11
Inputs types: [TensorType(int64, (True,)), TensorType(int32, vector), TensorType(float64, (True,)), TensorType(float64, (True,))]
Inputs shapes: [(1,), (1100,), (1,), (1,)]
Inputs strides: [(4,), (4,), (8,), (8,)]
Inputs values: [array([549]), 'not shown', array([ 1.07762995]), array([ 1.01502801])]
Outputs clients: [[Elemwise{eq,no_inplace}(Elemwise{Composite{Switch(GE(i0, i1), i2, i3)}}.0, TensorConstant{(1,) of 0}), Elemwise{Composite{Switch(GE(i0, i1), ((Switch(i2, i3, (i4 * log(i0))) - i5) - i0), i3)}}[(0, 0)](Elemwise{Composite{Switch(GE(i0, i1), i2, i3)}}.0, TensorConstant{(1,) of 0}, InplaceDimShuffle{x}.0, TensorConstant{(1,) of -inf}, TensorConstant{[ 13.  13...  0.   1.]}, TensorConstant{[ 22.55216...        ]})]]

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.

Being simple analyst, should i learn all this stuff about theano to being able to work with my statistical problems? Is a new mcmc sampler with gradient feature is only one thing that should motivates me to switch from pymc2 to pymc3?

Answer 1

For your first question, it looks like you're trying to pass a theano function as a variable. You need to call the function with the other variables as arguments, which will then return a theano variable. Try changing your line to

disasters=pymc3.Poisson('disasters', mu=rate1(switchpoint, mu1, mu2), observed = data)

I couldn't reproduce the error in your second part; the sampling worked just fine for me.