title: Using Categorical with multi-dimensional p in PyMC3
body: <p>I am running into problems when I am trying to use pm.Categorical to sample many instances at once (either with multidimensional p or using theano.scan). What is the best way to go here? My goal is to sample one response per draw for each of many individuals, based on the probability of each response for each individual.</p>
  
  <p>Here is a working example for many individuals, one draw:</p>
  
  <pre><code>import pymc3 as pm
  import theano.tensor as T
  import numpy as np
  
  actions = T.as_tensor_variable(np.array([0, 1, 2]))  # indiv0 selects response0, indiv1 response1, etc.
  with pm.Model() as model:
      p = pm.Beta('p', alpha=2, beta=2, shape=[3, 3])  # prob. for each indiv. for each response
      actions = pm.Categorical('actions', p=p, observed=actions)
      trace = pm.sample()
  </code></pre>
  
  <p>Everything works fine here. The model samples nine p's, one for each individual and each response. But what I'd like to do is to sample several responses for each individual. I first tried to add a dimension of trials to the actions data:</p>
  
  <pre><code>n_trials = 10
  actions = T.as_tensor_variable(np.array([[0, 1, 2]] * n_trials))  # same thing, 10 trials
  with pm.Model() as model:
      p = pm.Beta('p', alpha=2, beta=2, shape=[3, 3])  # same as above, but more trials
      actions = pm.Categorical('actions', p=p, observed=actions)
      trace = pm.sample()
  </code></pre>
  
  <p>Obviously, this leads to an error:</p>
  
  <pre><code>Wrong number of dimensions: expected 1, got 2 with shape (10, 3).
  </code></pre>
  
  <p>Surprisingly, it also leads to an error when I draw a separate p for each sample, each individual, and each response:</p>
  
  <pre><code>p = pm.Beta('p', alpha=2, beta=2, shape=[n_trials, 3, 3])
  </code></pre>
  
  <p>It seems like pm.Categorical only accepts 2-dimensional p. To keep p 2-dimensional, I tried scanning:</p>
  
  <pre><code>actions = T.as_tensor_variable(np.array([[0, 1, 2]] * 10)
  with pm.Model() as model:
      p = pm.Beta('p', alpha=2, beta=2, shape=[3, 3])  # same as above
      actions, _ = theano.scan(fn=lambda action, p: pm.Categorical('actions', p=p, observed=action),
                               sequences=[actions],
                               non_sequences=[p])
      trace = pm.sample()
  </code></pre>
  
  <p>The code fails with the following error message:</p>
  
  <pre><code>theano.gof.fg.MissingInputError: Input 0 of the graph (indices start from 0), used to compute Elemwise{Cast{int64}}(&lt;TensorType(int32, vector)&gt;), was not provided and not given a value. Use the Theano flag exception_verbosity='high', for more information on this error.
  </code></pre>
  
  <p>So maybe it is not possible to use pm.Categorical within theano.scan? Any help as to how to approach this problem would be super helpful!</p>