title: Softmax choice probabilities with Categorical in PyMC3
body: <p>I am trying to perform a parameter estimation for a single parameter of a softmax choice function in the following scenario:</p>
  
  <p>In each trial, three option values are given (e.g., [1,2,3]), and a subject makes a choice between the options (0, 1 or 2). The softmax function transforms option values into choice probabilities (vector of 3 probabilities, summing to 1), depending on a temperature parameter (here bound between 0 and 10).</p>
  
  <p>The choice in each trial is supposed to be modelled as a Categorical distribution with trial choice probabilities calculated from the softmax. Note that the choice probabilities of the Categorical depend on the option values and are therefore different in each trial.</p>
  
  <p>Here's what I came up with:</p>
  
  <pre class="lang-python prettyprint-override"><code># Generate data
  nTrials = 60 # number of trials (value triplets and choices)
  np.random.seed(42)
  # generate nTrials triplets of values
  values = np.random.choice([1,2,3,4,5], size=(nTrials, 3)) 
  
  choices = values.argmax(axis=1) # choose highest value option
  # add some random variation, so that *not* always the highest value option is chosen
  errors = np.random.rand(nTrials)&gt;0.8 # determine trials with non-optimal choice
  # randomly determine new choices for these trials
  choices[errors] = np.random.choice([0,1,2], size=sum(errors==True))
  
  # Model specification &amp; estimation
  import pymc3 as pm
  from theano import tensor as t
  with pm.Model():
  
      # prior over theta
      theta = pm.Uniform('theta', lower=0, upper=10)
  
      # softmax implementation
      enumerator  = pm.exp(theta*values) 
      denominator = t.reshape(pm.sum(pm.exp(theta*values), axis=1), (nTrials, 1))
      ps = enumerator/denominator
  
      # Likelihood (sampling model for the data)
      for trial in range(nTrials):
          yobs = pm.Categorical('yobs{}'.format(trial), p=ps[trial], observed=choices[trial])
  
      # draw 500 samples from posterior
      trace = pm.sample(500, pm.Metropolis())
  </code></pre>
  
  <p>This code fails for nTrials bigger than something like 50 with an extremely long warning / error message:</p>
  
  <p>Warning:</p>
  
  <pre><code>INFO (theano.gof.compilelock): Refreshing lock /Users/felixmolter/.theano/compiledir_Darwin-14.4.0-x86_64-i386-64bit-i386-2.7.8-64/lock_dir/lock
  INFO:theano.gof.compilelock:Refreshing lock /Users/felixmolter/.theano/compiledir_Darwin-14.4.0-x86_64-i386-64bit-i386-2.7.8-64/lock_dir/lock
  00001   #include &lt;Python.h&gt;
  00002   #include &lt;iostream&gt;
  00003   #include &lt;math.h&gt;
  00004   #include &lt;numpy/arrayobject.h&gt;
  00005   #include &lt;numpy/arrayscalars.h&gt;
  00006   #include &lt;vector&gt;
  00007   #include &lt;algorithm&gt;
  00008   //////////////////////
  00009   ////  Support Code
  00010   //////////////////////
  00011   
  00012   
  00013       namespace {
  00014       struct __struct_compiled_op_65734e56ae54d89bdcf84e36893358e6 {
  00015           PyObject* __ERROR;
  00016   
  00017           PyObject* storage_V3;
  00018   PyObject* storage_V5;
  00019   PyObject* storage_V7;
  00020   PyObject* storage_V9;
  00021   PyObject* storage_V11;
  00022   PyObject* storage_V13;
  [...]
  </code></pre>
  
  <p>Error:</p>
  
  <pre><code>Exception: ('The following error happened while compiling the node', Elemwise{Composite{((Switch(LE(Abs((i0 + i1)), i2), log(i3), i4) + Switch(LE(Abs((i0 + i5)), i2), log(i6), i4) + Switch(LE(Abs((i0 + i7)), i2), log(i8), i4) + Switch(LE(Abs((i0 + i9)), i2), log(i10), i4) + Switch(LE(Abs((i0 + i11)), [...]
  </code></pre>
  
  <p>I am rather new to PyMC (and Theano) and I feel my implementation is really clunky and suboptimal. Any help and advice is strongly appreciated!</p>
  
  <p>Felix</p>
  
  <p>Edit: I've uploaded the code as a notebook, showing the warnings and error messages in full: <a href="http://nbviewer.ipython.org/github/moltaire/softmaxPyMC/blob/master/softmax_stackoverflow.ipynb" rel="nofollow">http://nbviewer.ipython.org/github/moltaire/softmaxPyMC/blob/master/softmax_stackoverflow.ipynb</a></p>