stick breaking transform does not work for 2D vars

[twiecki]

from pymc3 import Dirichlet, Model
with Model():
    Dirichlet('beta', np.ones((3, 3)))

produces:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-15-059ef57ab101> in <module>()
      1 from pymc3 import Dirichlet, Model
      2 with Model():
----> 3     Dirichlet('beta', np.ones((3, 3)))

/home/wiecki/working/projects/pymc/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     17             data = kwargs.pop('observed', None)
     18             dist = cls.dist(*args, **kwargs)
---> 19             return model.Var(name, dist, data)
     20         elif name is None:
     21             return object.__new__(cls) #for pickle

/home/wiecki/working/projects/pymc/pymc3/model.py in Var(self, name, dist, data)
    151                 self.free_RVs.append(var)
    152             else:
--> 153                 var = TransformedRV(name=name, distribution=dist, model=self, transform=dist.transform)
    154                 self.deterministics.append(var)
    155                 return var

/home/wiecki/working/projects/pymc/pymc3/model.py in __init__(self, type, owner, index, name, distribution, model, transform)
    506             self.model = model
    507 
--> 508             self.transformed = model.Var(name + "_" + transform.name, transform.apply(distribution))
    509 
    510             normalRV = transform.backward(self.transformed)

/home/wiecki/working/projects/pymc/pymc3/distributions/transforms.py in apply(self, dist)
     23 
     24     def apply(self, dist):
---> 25         return TransformedDistribution.dist(dist, self)
     26 
     27     def __str__(self):

/home/wiecki/working/projects/pymc/pymc3/distributions/distribution.py in dist(cls, *args, **kwargs)
     29     def dist(cls, *args, **kwargs):
     30         dist = object.__new__(cls)
---> 31         dist.__init__(*args, **kwargs)
     32         return dist
     33 

/home/wiecki/working/projects/pymc/pymc3/distributions/transforms.py in __init__(self, dist, transform, *args, **kwargs)
     46             arguments to Distribution"""
     47         forward = transform.forward
---> 48         testval = forward(dist.default())
     49 
     50         self.dist = dist

/home/wiecki/working/projects/pymc/pymc3/distributions/transforms.py in forward(self, x)
    151         #reverse cumsum
    152         x0 = x[:-1]
--> 153         s = t.extra_ops.cumsum(x0[::-1])[::-1] + x[-1]
    154         z = x0/s
    155         Km1 = x.shape[0] - 1

/home/wiecki/envs/pymc3/local/lib/python2.7/site-packages/theano/tensor/var.py in __add__(self, other)
    126     def __add__(self, other):
    127         try:
--> 128             return theano.tensor.basic.add(self, other)
    129         # We should catch the minimum number of exception here.
    130         # Otherwise this will convert error when Theano flags

/home/wiecki/envs/pymc3/local/lib/python2.7/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    556                 thunk.outputs = [storage_map[v] for v in node.outputs]
    557 
--> 558                 required = thunk()
    559                 assert not required  # We provided all inputs
    560 

/home/wiecki/envs/pymc3/local/lib/python2.7/site-packages/theano/gof/op.py in rval()
    733 
    734         def rval():
--> 735             fill_storage()
    736             for o in node.outputs:
    737                 compute_map[o][0] = True

/home/wiecki/envs/pymc3/local/lib/python2.7/site-packages/theano/gof/cc.py in __call__(self)
   1538                 raise
   1539 
-> 1540             raise exc_type, exc_value, exc_trace
   1541 
   1542 

ValueError: Input dimension mis-match. (input[0].shape[0] = 6, input[1].shape[0] = 3)

CC @jsalvatier

[twiecki]

I think we need to rewrite forward() to work on the last axis:

    def forward(self, x):
        #reverse cumsum
        x0 = x[..., :-1]
        s = t.extra_ops.cumsum(x0[::-1], axis=-1)[::-1] + x[...,[-1]]
        z = x0/s
        Km1 = x.shape[-1] - 1
        k = arange(Km1)
        eq_share = - t.log(Km1 - k) # logit(1./(Km1 + 1 - k))
        y = logit(z) - eq_share
        return y

This doesn't quite work yet though.

[jsalvatier]

Good point. Yeah, this fits in with multivariate distributions needing a special axis.

[fonnesbeck]

@twiecki is this fixed by #844 ?

[twiecki]

Yes, I believe so.