Issue with shape of GaussianRandomWalk variable

Questions
1

I’m interested in running a simple time series model with the GaussianRandomWalk distribution. My problem is that I appear to be unable to specify the shape of a GRW variable correctly. In the code below, ‘observed’ is a vector with shape (31,).

    t = len(observed)
    with pymc3.Model() as random_walk_model:
        alpha = pymc3.InverseGamma('alpha',alpha = 1.0,beta = 1.0)
        sd    = pymc3.InverseGamma('sd',alpha = 2.0,beta = 2.0)
        s     = pymc3.GaussianRandomWalk('s',sd = sd,shape = t)
        x     = pymc3.Weibull('x',alpha = alpha,beta = s,shape = t,observed = observed)
        trace = pymc3.sample()

Upon attempting to initialize the variable ‘x’, I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-560-0377f05d8ad5> in <module>()
      6     s     = pymc3.GaussianRandomWalk('s',sd = sd,shape = t)
      7     print s.dshape
----> 8     x     = pymc3.Weibull('x',alpha = alpha,beta = s,shape = t,observed = observed)
      9     trace = pymc3.sample() 
     10 

/home/ubuntu/anaconda2/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
     34                 raise TypeError("observed needs to be data but got: {}".format(type(data)))
     35             total_size = kwargs.pop('total_size', None)
---> 36             dist = cls.dist(*args, **kwargs)
     37             return model.Var(name, dist, data, total_size)
     38         else:

/home/ubuntu/anaconda2/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in dist(cls, *args, **kwargs)
     45     def dist(cls, *args, **kwargs):
     46         dist = object.__new__(cls)
---> 47         dist.__init__(*args, **kwargs)
     48         return dist
     49 

/home/ubuntu/anaconda2/lib/python2.7/site-packages/pymc3/distributions/continuous.pyc in __init__(self, alpha, beta, *args, **kwargs)
   1499 
   1500         assert_negative_support(alpha, 'alpha', 'Weibull')
-> 1501         assert_negative_support(beta, 'beta', 'Weibull')
   1502 
   1503     def random(self, point=None, size=None, repeat=None):

/home/ubuntu/anaconda2/lib/python2.7/site-packages/pymc3/distributions/continuous.pyc in assert_negative_support(var, label, distname, value)
     58         try:
     59             # Untransformed distribution
---> 60             support = np.isfinite(var.distribution.logp(value).tag.test_value)
     61         except AttributeError:
     62             # Otherwise no direct evidence of non-positive support

/home/ubuntu/anaconda2/lib/python2.7/site-packages/pymc3/distributions/timeseries.pyc in logp(self, x)
    161         init = self.init
    162 
--> 163         x_im1 = x[:-1]
    164         x_i = x[1:]
    165 

TypeError: 'float' object has no attribute '__getitem__'

I’ve tried playing around with the dimension of the variables here but can’t seem to get it to work. Any advice here?

2

When you are defining an observed variable, you don’t need to specify the shape kwarg.

3

Thanks, that’s good to know - doesn’t fix the resulting error though.

4

Right, I think the Weibull distribution does not support getting theano tensor as input, which is weird.

For now, you can cherry pick the logp from pm.Weibull and define a pm.DensityDist:

import pymc3 as pm
import numpy as np
import scipy.stats as st
import theano.tensor as tt
bound=pm.distributions.dist_math.bound
alpha, beta = 1.5, 2.
observed = st.weibull_min(alpha, scale=beta).rvs(30)
t = len(observed)
def Weibull_logp(alpha, beta, value):
    return bound(tt.log(alpha) - tt.log(beta)
                  + (alpha - 1) * tt.log(value / beta)
                  - (value / beta)**alpha,
                     value >= 0, alpha > 0, beta > 0)
with pm.Model():
    alpha = pm.InverseGamma('alpha',alpha = 1.0,beta = 1.0)
    sd    = pm.InverseGamma('sd',alpha = 2.0,beta = 2.0)
    s     = pm.GaussianRandomWalk('s',sd = sd,shape = t)
    x     = pm.DensityDist('x', Weibull_logp, observed=dict(alpha=alpha, beta=sd, value=observed))
    trace = pm.sample()
5

I’ll try this out. Thanks a ton!

1 Like