pm.Categorical behaves differently in a model versus as pm.Categorical.dist

[MariaEckstein]

Using pm.Categorical.dist using a 2-d array of probabilities samples multiple values:

p = 1 / 3 * np.ones([30, 3])
pm.Categorical.dist(p=p).random(1)
>>> array([2, 1, 0, 0, 1, 2, 2, 2, 0, 0, 0, 1, 0, 0, 0, 2, 0, 1, 0, 0, 1, 0,
       0, 2, 1, 1, 2, 1, 0, 1], dtype=int64)

The same thing within a model samples just one value:

p = 1 / 3 * np.ones([30, 3])
with pm.Model() as model:
    p = T.as_tensor_variable(p)
    choice = pm.Categorical('choice', p=p)
    T.printing.Print('choice')(choice)  # prints choice __str__ = 0
    trace = pm.sample(100, tune=10, chains=1, cores=1)

Interestingly, when passed an array of length 30 as observed, pm.Categorical draws multiple values again:

with pm.Model() as model:
    p = pm.Uniform('p', lower=0, upper=1, shape=3)
    p_tile = T.tile(p, 30).reshape([30, 3])
    choice = pm.Categorical('choice', p=p_tile, observed=np.array([0, 0, 2] * 10))
    T.printing.Print('choice')(choice)  # prints choice __str__ = [0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2]
    trace = pm.sample(5000, tune=500, chains=2, cores=1)
print(pm.summary(trace))

>>>
                     mean        sd  mc_error   hpd_2.5  hpd_97.5
p__0             0.746065  0.196305  0.013013  0.321824  0.995407  # estimates of p correspond to observed
p__1             0.043839  0.051256  0.002698  0.000105  0.140824
p__2             0.460127  0.208897  0.012058  0.107075  0.903796

The same model without observations seems to be doing something quite different:

with pm.Model() as model:
    p = pm.Uniform('p', lower=0, upper=1, shape=3)
    p_tile = T.tile(p, 30).reshape([30, 3])
    choice = pm.Categorical('choice', p=p_tile)
    T.printing.Print('choice')(choice)  # prints choice __str__ = 0
    trace = pm.sample(500, tune=500, chains=1, cores=1)
print(pm.summary(trace))

>>>
                 mean        sd  mc_error   hpd_2.5  hpd_97.5
choice       1.000000  0.000000  0.000000  1.000000  1.000000  # don't understand where these estimates 
p__0         0.026698  0.031869  0.001404  0.000109  0.085003
p__1         0.743670  0.184351  0.009447  0.377753  0.999541
p__2         0.026075  0.028349  0.001499  0.000077  0.081620

• PyMC3 Version: 3.4.1
• Theano Version: 1.0.2
• Python Version: 3.6
• Operating system: Windows

[junpenglao]

Shape of the RV is not broadcast from parameters for Categorical and Multinomial. Note that in both case if you pass the shape=... to the RV the result would be correct.

In [5]: p = 1 / 3 * np.ones([30, 3])
   ...: with pm.Model() as m:
   ...:     choice = pm.Categorical('choice', p=p)
   ...:     trace = pm.sample(100, tune=0, chains=1)
   ...: trace['choice'].shape # should be (100, 30)

Only 100 samples in chain.
Sequential sampling (1 chains in 1 job)
CategoricalGibbsMetropolis: [choice]
100%|███████████████████████████████████████| 100/100 [00:00<00:00, 5988.87it/s]
Only one chain was sampled, this makes it impossible to run some convergence checks
Out[5]: (100,)

In [6]: p = 1 / 3 * np.ones([30, 3])
   ...: with pm.Model() as m:
   ...:     choice = pm.Multinomial('choice', n=1, p=p)
   ...:     trace = pm.sample(100, tune=0, chains=1)
   ...: trace['choice'].shape # should be (100, 30)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-6-5177f41cfeba> in <module>()
      1 p = 1 / 3 * np.ones([30, 3])
      2 with pm.Model() as m:
----> 3     choice = pm.Multinomial('choice', n=1, p=p)
      4     trace = pm.sample(100, tune=0, chains=1)
      5 trace['choice'].shape # should be (100, 30)

~/Documents/pymc3/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     40             total_size = kwargs.pop('total_size', None)
     41             dist = cls.dist(*args, **kwargs)
---> 42             return model.Var(name, dist, data, total_size)
     43         else:
     44             raise TypeError("Name needs to be a string but got: {}".format(name))

~/Documents/pymc3/pymc3/model.py in Var(self, name, dist, data, total_size)
    806                 with self:
    807                     var = FreeRV(name=name, distribution=dist,
--> 808                                  total_size=total_size, model=self)
    809                 self.free_RVs.append(var)
    810             else:

~/Documents/pymc3/pymc3/model.py in __init__(self, type, owner, index, name, distribution, total_size, model)
   1206             self.tag.test_value = np.ones(
   1207                 distribution.shape, distribution.dtype) * distribution.default()
-> 1208             self.logp_elemwiset = distribution.logp(self)
   1209             # The logp might need scaling in minibatches.
   1210             # This is done in `Factor`.

~/Documents/pymc3/pymc3/distributions/multivariate.py in logp(self, x)
    579 
    580         return bound(
--> 581             factln(n) + tt.sum(-factln(x) + logpow(p, x), axis=-1, keepdims=True),
    582             tt.all(x >= 0),
    583             tt.all(tt.eq(tt.sum(x, axis=-1, keepdims=True), n)),

~/Documents/pymc3/pymc3/distributions/dist_math.py in factln(n)
     70 
     71 def factln(n):
---> 72     return gammaln(n + 1)
     73 
     74 

~/anaconda3/lib/python3.6/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

~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    623             for i, ins in enumerate(node.inputs):
    624                 try:
--> 625                     storage_map[ins] = [self._get_test_value(ins)]
    626                     compute_map[ins] = [True]
    627                 except AttributeError:

~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in _get_test_value(cls, v)
    560             # ensure that the test value is correct
    561             try:
--> 562                 ret = v.type.filter(v.tag.test_value)
    563             except Exception as e:
    564                 # Better error message.

~/anaconda3/lib/python3.6/site-packages/theano/tensor/type.py in filter(self, data, strict, allow_downcast)
    176             raise TypeError("Wrong number of dimensions: expected %s,"
    177                             " got %s with shape %s." % (self.ndim, data.ndim,
--> 178                                                         data.shape))
    179         if not data.flags.aligned:
    180             try:

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 2 with shape (30, 3).

[lucianopaz]

#3386 will make the case where no shape is provided to the Categorical RV to fail with the following TypeError:

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 1 with shape (30,).

Furthermore, it will also make the following snippet give a bad initial energy, because p will not be normalized by default when it's a TensorVariable

with pm.Model() as model:
    p = pm.Uniform('p', lower=0, upper=1, shape=3)
    p_tile = T.tile(p, 30).reshape([30, 3])
    choice = pm.Categorical('choice', p=p_tile, observed=np.array([0, 0, 2] * 10))
    T.printing.Print('choice')(choice)  # prints choice __str__ = [0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2 0 0 2]
    trace = pm.sample(5000, tune=500, chains=2, cores=1)
print(pm.summary(trace))

This means that #3386, will remove the ambiguity, but this issue also highlights the potential backwards compatibility dangers we might run into if the normalization is not applied in the logp

[twiecki]

Closed by #3386.