Both the logistic regression analyses of subjects choices, and reinforcement learning model fitting used a Bayesian hierarchical modeling framework (Huys et al., 2011), in which parameter vectors hi for individual sessions were assumed to be drawn from Gaussian distributions at the population level with means and variance θ={μ,Σ}. The population level prior distributions were set to their maximum likelihood estimate:19θML=argmaxθ{p(D|θ)=argmaxθ{∏iN∫dhip(Di|hi)p(hi|θ)}Optimization was performed using the Expectation-Maximization algorithm with a Laplace approximation for the E-step at the k-th iteration given by:20p(hik|Di)=N(mik,Vik)21mik=argmaxh{p(Di|h)p(h|θk−1)}Where N(mik,Vik) is a normal distribution with mean mik given by the maximum a posteriori value of the session parameter vector hi given the population level means and variance θk−1, and the covariance Vik given by the inverse Hessian of the likelihood around mik. For simplicity we assumed that the population level covariance Σ had zero off-diagonal terms. For the k-th M-step of the EM algorithm the population level prior distribution parameters θ={μ,Σ} are updated as:22μk=1N∑i=1Nmik23Σ=1N∑i=1N[(mik)2+Vik]−(μk)2Parameters were transformed before inference to enforce constraints (0<{Gmf,Gmo,Gmb},0<{αQ,fQ,λ,αT,fT,αc,αm}<1).
