In [1]:
import theano
theano.config.exception_verbosity='high'
theano.config.mode='DebugMode'
import pymc3 as pm

%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
import seaborn as sns
sns.set_style("ticks")
sns.set_context("poster")

In [2]:
disaster_data = np.ma.masked_values(np.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, -999, 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, -999, 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], dtype=np.int64), value=-999,)
year = np.arange(1851, 1962)

In [3]:
fig, ax = plt.subplots(1,2)

ax[0].hist(disaster_data, bins=range(disaster_data.max() + 1))
ax[1].plot(disaster_data, "o")
ax[1].set_xlabel("Year")
ax[1].set_ylabel("Disaster Count")
ax[0].set_ylabel("Disaster Count")


Out[3]:
<matplotlib.text.Text at 0xd19aa20>

In [4]:
disaster_data


Out[4]:
masked_array(data = [4L 5L 4L 0L 1L 4L 3L 4L 0L 6L 3L 3L 4L 0L 2L 6L 3L 3L 5L 4L 5L 3L 1L 4L 4L
 1L 5L 5L 3L 4L 2L 5L 2L 2L 3L 4L 2L 1L 3L -- 2L 1L 1L 1L 1L 3L 0L 0L 1L 0L
 1L 1L 0L 0L 3L 1L 0L 3L 2L 2L 0L 1L 1L 1L 0L 1L 0L 1L 0L 0L 0L 2L 1L 0L 0L
 0L 1L 1L 0L 2L 3L 3L 1L -- 2L 1L 1L 1L 1L 2L 4L 2L 0L 0L 1L 4L 0L 0L 0L 1L
 0L 0L 0L 0L 0L 1L 0L 0L 1L 0L 1L],
             mask = [False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False  True False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False False False False False False False False False  True
 False False False False False False False False False False False False
 False False False False False False False False False False False False
 False False False],
       fill_value = -999)

In [5]:
year


Out[5]:
array([1851, 1852, 1853, 1854, 1855, 1856, 1857, 1858, 1859, 1860, 1861,
       1862, 1863, 1864, 1865, 1866, 1867, 1868, 1869, 1870, 1871, 1872,
       1873, 1874, 1875, 1876, 1877, 1878, 1879, 1880, 1881, 1882, 1883,
       1884, 1885, 1886, 1887, 1888, 1889, 1890, 1891, 1892, 1893, 1894,
       1895, 1896, 1897, 1898, 1899, 1900, 1901, 1902, 1903, 1904, 1905,
       1906, 1907, 1908, 1909, 1910, 1911, 1912, 1913, 1914, 1915, 1916,
       1917, 1918, 1919, 1920, 1921, 1922, 1923, 1924, 1925, 1926, 1927,
       1928, 1929, 1930, 1931, 1932, 1933, 1934, 1935, 1936, 1937, 1938,
       1939, 1940, 1941, 1942, 1943, 1944, 1945, 1946, 1947, 1948, 1949,
       1950, 1951, 1952, 1953, 1954, 1955, 1956, 1957, 1958, 1959, 1960,
       1961])

In [7]:
with pm.Model() as disaster_model:

    switchpoint = pm.DiscreteUniform('switchpoint', lower=year.min(), upper=year.max(), testval=1900)

    # Priors for pre- and post-switch rates number of disasters
    early_rate = pm.Exponential('early_rate', 1)
    late_rate = pm.Exponential('late_rate', 1)

    # Allocate appropriate Poisson rates to years before and after current
    rate = pm.math.switch(switchpoint >= year, early_rate, late_rate)
    

    disasters = pm.Poisson('disasters', rate, observed=disaster_data)

In [8]:
with disaster_model:
    step1 = pm.NUTS([early_rate, late_rate])

    # Use Metropolis for switchpoint, and missing values since it accommodates discrete variables
    step2 = pm.Metropolis([switchpoint, disasters.missing_values[0]] )

    trace = pm.sample(10000, step=[step1, step2])


---------------------------------------------------------------------------
InvalidValueError                         Traceback (most recent call last)
<ipython-input-8-3a4f2e83047f> in <module>()
      1 with disaster_model:
----> 2     step1 = pm.NUTS([early_rate, late_rate])
      3 
      4     # Use Metropolis for switchpoint, and missing values since it accommodates discrete variables
      5     step2 = pm.Metropolis([switchpoint, disasters.missing_values[0]] )

E:\Anaconda\lib\site-packages\pymc3\step_methods\nuts.pyc in __init__(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)
     70         if isinstance(scaling, dict):
     71             scaling = guess_scaling(
---> 72                 Point(scaling, model=model), model=model, vars=vars)
     73 
     74         n = scaling.shape[0]

E:\Anaconda\lib\site-packages\pymc3\tuning\scaling.pyc in guess_scaling(point, vars, model, scaling_bound)
    109     model = modelcontext(model)
    110     try:
--> 111         h = find_hessian_diag(point, vars, model=model)
    112     except NotImplementedError:
    113         h = fixed_hessian(point, vars, model=model)

E:\Anaconda\lib\site-packages\pymc3\tuning\scaling.pyc in find_hessian_diag(point, vars, model)
    102     """
    103     model = modelcontext(model)
--> 104     H = model.fastfn(hessian_diag(model.logpt, vars))
    105     return H(Point(point, model=model))
    106 

E:\Anaconda\lib\site-packages\pymc3\memoize.pyc in memoizer(*args, **kwargs)
     13 
     14         if key not in cache:
---> 15             cache[key] = obj(*args, **kwargs)
     16 
     17         return cache[key]

E:\Anaconda\lib\site-packages\pymc3\theanof.pyc in hessian_diag(f, vars)
    107 
    108     if vars:
--> 109         return -tt.concatenate([hessian_diag1(f, v) for v in vars], axis=0)
    110     else:
    111         return empty_gradient

E:\Anaconda\lib\site-packages\pymc3\theanof.pyc in hessian_diag1(f, v)
     98         return gradient1(g[i], v)[i]
     99 
--> 100     return theano.map(hess_ii, idx)[0]
    101 
    102 

E:\Anaconda\lib\site-packages\theano\scan_module\scan_views.pyc in map(fn, sequences, non_sequences, truncate_gradient, go_backwards, mode, name)
     68                      go_backwards=go_backwards,
     69                      mode=mode,
---> 70                      name=name)
     71 
     72 

E:\Anaconda\lib\site-packages\theano\scan_module\scan.pyc in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict)
   1064             pass
   1065         scan_inputs += [arg]
-> 1066     scan_outs = local_op(*scan_inputs)
   1067     if type(scan_outs) not in (list, tuple):
   1068         scan_outs = [scan_outs]

E:\Anaconda\lib\site-packages\theano\gof\op.pyc in __call__(self, *inputs, **kwargs)
    661                 thunk.outputs = [storage_map[v] for v in node.outputs]
    662 
--> 663                 required = thunk()
    664                 assert not required  # We provided all inputs
    665 

E:\Anaconda\lib\site-packages\theano\scan_module\scan_op.pyc in rval(p, i, o, n, allow_gc)
    951         def rval(p=p, i=node_input_storage, o=node_output_storage, n=node,
    952                  allow_gc=allow_gc):
--> 953             r = p(n, [x[0] for x in i], o)
    954             for o in node.outputs:
    955                 compute_map[o][0] = True

E:\Anaconda\lib\site-packages\theano\scan_module\scan_op.pyc in p(node, args, outs)
    940                                                 args,
    941                                                 outs,
--> 942                                                 self, node)
    943         except (ImportError, theano.gof.cmodule.MissingGXX):
    944             p = self.execute

theano/scan_module/scan_perform.pyx in theano.scan_module.scan_perform.perform (C:\Users\Shubhanshu\AppData\Local\Theano\compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_58_Stepping_9_GenuineIntel-2.7.12-64\scan_perform\mod.cpp:4415)()

theano/scan_module/scan_perform.pyx in theano.scan_module.scan_perform.perform (C:\Users\Shubhanshu\AppData\Local\Theano\compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_58_Stepping_9_GenuineIntel-2.7.12-64\scan_perform\mod.cpp:4199)()

E:\Anaconda\lib\site-packages\theano\compile\debugmode.pyc in deco()
   2303                     self.maker.mode.check_isfinite
   2304                 try:
-> 2305                     return f()
   2306                 finally:
   2307                     # put back the filter_checks_isfinite

E:\Anaconda\lib\site-packages\theano\compile\debugmode.pyc in f()
   1960                                 r, storage_map[r][0],
   1961                                 hint=("Graph Input '%s' has invalid value "
-> 1962                                       "%s" % (r, storage_map[r][0])))
   1963                         r_vals[r] = storage_map[r][0]
   1964                         storage_map[r][0] = None

InvalidValueError: InvalidValueError
        type(variable) = TensorType(float32, scalar)
        variable       = TensorConstant{-inf}
        type(value)    = <type 'numpy.ndarray'>
        dtype(value)   = float32
        shape(value)   = ()
        value          = -inf
        min(value)     = -inf
        max(value)     = -inf
        isfinite       = False
        client_node    = None
        hint           = Graph Input 'TensorConstant{-inf}' has invalid value -inf
        specific_hint  = none
        context        = ...
  TensorConstant{-inf} [id A]

        

In [9]:
with disaster_model:
    trace = pm.sample(100)


Assigned Metropolis to switchpoint
Assigned NUTS to early_rate_log_
Assigned NUTS to late_rate_log_
Assigned Metropolis to disasters_missing
---------------------------------------------------------------------------
InvalidValueError                         Traceback (most recent call last)
<ipython-input-9-92b5c49166ad> in <module>()
      1 with disaster_model:
----> 2     trace = pm.sample(100)

E:\Anaconda\lib\site-packages\pymc3\sampling.pyc in sample(draws, step, init, n_init, start, trace, chain, njobs, tune, progressbar, model, random_seed)
    151             start = start_
    152     else:
--> 153         step = assign_step_methods(model, step)
    154 
    155     if njobs is None:

E:\Anaconda\lib\site-packages\pymc3\sampling.pyc in assign_step_methods(model, step, methods)
     74 
     75     # Instantiate all selected step methods
---> 76     steps += [step(vars=selected_steps[step]) for step in selected_steps if selected_steps[step]]
     77 
     78     if len(steps) == 1:

E:\Anaconda\lib\site-packages\pymc3\step_methods\nuts.pyc in __init__(self, vars, scaling, step_scale, is_cov, state, Emax, target_accept, gamma, k, t0, model, profile, **kwargs)
     70         if isinstance(scaling, dict):
     71             scaling = guess_scaling(
---> 72                 Point(scaling, model=model), model=model, vars=vars)
     73 
     74         n = scaling.shape[0]

E:\Anaconda\lib\site-packages\pymc3\tuning\scaling.pyc in guess_scaling(point, vars, model, scaling_bound)
    109     model = modelcontext(model)
    110     try:
--> 111         h = find_hessian_diag(point, vars, model=model)
    112     except NotImplementedError:
    113         h = fixed_hessian(point, vars, model=model)

E:\Anaconda\lib\site-packages\pymc3\tuning\scaling.pyc in find_hessian_diag(point, vars, model)
    102     """
    103     model = modelcontext(model)
--> 104     H = model.fastfn(hessian_diag(model.logpt, vars))
    105     return H(Point(point, model=model))
    106 

E:\Anaconda\lib\site-packages\pymc3\memoize.pyc in memoizer(*args, **kwargs)
     13 
     14         if key not in cache:
---> 15             cache[key] = obj(*args, **kwargs)
     16 
     17         return cache[key]

E:\Anaconda\lib\site-packages\pymc3\theanof.pyc in hessian_diag(f, vars)
    107 
    108     if vars:
--> 109         return -tt.concatenate([hessian_diag1(f, v) for v in vars], axis=0)
    110     else:
    111         return empty_gradient

E:\Anaconda\lib\site-packages\pymc3\theanof.pyc in hessian_diag1(f, v)
     98         return gradient1(g[i], v)[i]
     99 
--> 100     return theano.map(hess_ii, idx)[0]
    101 
    102 

E:\Anaconda\lib\site-packages\theano\scan_module\scan_views.pyc in map(fn, sequences, non_sequences, truncate_gradient, go_backwards, mode, name)
     68                      go_backwards=go_backwards,
     69                      mode=mode,
---> 70                      name=name)
     71 
     72 

E:\Anaconda\lib\site-packages\theano\scan_module\scan.pyc in scan(fn, sequences, outputs_info, non_sequences, n_steps, truncate_gradient, go_backwards, mode, name, profile, allow_gc, strict)
   1064             pass
   1065         scan_inputs += [arg]
-> 1066     scan_outs = local_op(*scan_inputs)
   1067     if type(scan_outs) not in (list, tuple):
   1068         scan_outs = [scan_outs]

E:\Anaconda\lib\site-packages\theano\gof\op.pyc in __call__(self, *inputs, **kwargs)
    661                 thunk.outputs = [storage_map[v] for v in node.outputs]
    662 
--> 663                 required = thunk()
    664                 assert not required  # We provided all inputs
    665 

E:\Anaconda\lib\site-packages\theano\scan_module\scan_op.pyc in rval(p, i, o, n, allow_gc)
    951         def rval(p=p, i=node_input_storage, o=node_output_storage, n=node,
    952                  allow_gc=allow_gc):
--> 953             r = p(n, [x[0] for x in i], o)
    954             for o in node.outputs:
    955                 compute_map[o][0] = True

E:\Anaconda\lib\site-packages\theano\scan_module\scan_op.pyc in p(node, args, outs)
    940                                                 args,
    941                                                 outs,
--> 942                                                 self, node)
    943         except (ImportError, theano.gof.cmodule.MissingGXX):
    944             p = self.execute

theano/scan_module/scan_perform.pyx in theano.scan_module.scan_perform.perform (C:\Users\Shubhanshu\AppData\Local\Theano\compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_58_Stepping_9_GenuineIntel-2.7.12-64\scan_perform\mod.cpp:4415)()

theano/scan_module/scan_perform.pyx in theano.scan_module.scan_perform.perform (C:\Users\Shubhanshu\AppData\Local\Theano\compiledir_Windows-10-10.0.14393-Intel64_Family_6_Model_58_Stepping_9_GenuineIntel-2.7.12-64\scan_perform\mod.cpp:4199)()

E:\Anaconda\lib\site-packages\theano\compile\debugmode.pyc in deco()
   2303                     self.maker.mode.check_isfinite
   2304                 try:
-> 2305                     return f()
   2306                 finally:
   2307                     # put back the filter_checks_isfinite

E:\Anaconda\lib\site-packages\theano\compile\debugmode.pyc in f()
   1960                                 r, storage_map[r][0],
   1961                                 hint=("Graph Input '%s' has invalid value "
-> 1962                                       "%s" % (r, storage_map[r][0])))
   1963                         r_vals[r] = storage_map[r][0]
   1964                         storage_map[r][0] = None

InvalidValueError: InvalidValueError
        type(variable) = TensorType(float32, scalar)
        variable       = TensorConstant{-inf}
        type(value)    = <type 'numpy.ndarray'>
        dtype(value)   = float32
        shape(value)   = ()
        value          = -inf
        min(value)     = -inf
        max(value)     = -inf
        isfinite       = False
        client_node    = None
        hint           = Graph Input 'TensorConstant{-inf}' has invalid value -inf
        specific_hint  = none
        context        = ...
  TensorConstant{-inf} [id A]

        

In [ ]: