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
