Softmax choice probabilities with Categorical in PyMC3
Question: I am trying to perform a parameter estimation for a single parameter of a
softmax choice function in the following scenario:

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).

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.

Here's what I came up with:

    
    
    # 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)>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 & 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())
    

This code fails for nTrials bigger than something like 50 with an extremely
long warning / error message:

Warning:

    
    
    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 <Python.h>
    00002   #include <iostream>
    00003   #include <math.h>
    00004   #include <numpy/arrayobject.h>
    00005   #include <numpy/arrayscalars.h>
    00006   #include <vector>
    00007   #include <algorithm>
    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;
    [...]
    

Error:

    
    
    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)), [...]
    

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!

Felix

Edit: I've uploaded the code as a notebook, showing the warnings and error
messages in full:
<http://nbviewer.ipython.org/github/moltaire/softmaxPyMC/blob/master/softmax_stackoverflow.ipynb>


Answer: I found again this case. Just as a follow up, now not able able to reproduce
it now. So I think it got fixed. We fixed related problem that could hve
caused in some case this error.

It work with g++ 4.5.1. If you have this problem update Theano to the
development version. If that don't fix it, try to use a g++ more recent, this
could be related to older g++ version.

