Many Divergences for a simple model

If I am reading your first post correctly, the observations are proportions (i.e., proportion of money given to person 1). In that case, I would probably use a Binomial distribution, which takes a probability (such as the sigmoid quantity returned from the logistic function) and a count (i.e., n), which would be the total number of units (dollars, quarters, dimes, pennies) that a user could give to person 1. Note that treating your data in this way requires that “unconvert” everything from the proportions you currently have to just the numerators of those proportions (and then n is the denominator). So instead of observing 30% of $10 dollars given to person 1, your observation would be that $3 out of a possible $10 were given to person 1 (or, equivalently, 300 cents out of a possible 1,000 cents).

Regarding the use of tt.repeat(), it may work, but it is a rather brittle way of handling subject-specific parameters. For example, if the number of trials per subject were to change or subjects’ data were intermixed in train_x/train_y (e.g., the first and third and fourth observations were from one subject and the second, fifth and sixth were from another), it wouldn’t be obvious that the call to tt.repeat() would need to be altered. I have mocked up a more conventional way of constructing and then referencing parameter arrays. Note that the subject_id list is zero-indexed (first subject is subject 0) and subjects are numbered consecutively. You can do this manually, but I usually convert subject identifiers into categorical variables via pandas and then use the category codes as my indicies into arrays of subject-specific parameters.

import numpy as np
import pymc3 as pm
import theano.tensor as tt

train_x = np.array([ 0.7837751 ,  0.19241506, -0.2124012 , -0.82463659,  0.78311322, 0.39562878, -1.13989023, -0.86822927, -0.60025648,  0.3742629 ])
train_y = np.array([5.63498758e-01, 2.88627774e-01, 2.92535202e-01, 3.35866828e-01, 6.67700455e-01, 9.99999856e-01, 1.45195098e-07, 1.34094597e-07, 1.02634229e-07, 9.99999885e-01])
hold_out_x = np.array([-1.47333049,  0.19393238, -0.21990887,  0.74694174, -0.34769498, 0.26646929,  0.72834725, -0.47401384, -1.11955062,  0.46028331])
hold_out_y =  np.array([0.21395213, 0.41452217, 0.47663821, 0.79562313, 0.08839646, 0.62057976, 0.765773  , 0.28458744, 0.1346239 , 0.69640536])

###########
# new stuff

# assume the first 5 observations are from subject 0
# and the next 5 are from subject 1
subject_id = [0,0,0,0,0,1,1,1,1,1]

# assume that the amount of money potentially given to person 1 is 1000
n_pennies = 1000
# you could instead allow different observations to have different amounts
# in that case you could do something like this:
# n_pennies = [1000, 500, 100, 2000, 300, 100, 250, 700, 350, 1200]

# now discretize the observations so that they are integers (e.g., how many pennies were given to person 1) 
train_y_discrete = (n_pennies * train_y).astype(int)

# end (most of the) new stuff
###########


shape = len(train_x)
person_shape = int(shape/5)

with pm.Model() as model_others:
    # this isn't necessary
    #slopes = pm.Data("slopes", train_x)
    mu = pm.Normal('mu', mu=0, sigma=1)
    sigma = pm.HalfNormal('sigma', sd=2.0) 
    beta = pm.Normal('beta', mu=mu, sigma=sigma, shape=person_shape)

    # here, we use subject_id to select the appropriate beta
    p = pm.Deterministic('p', 1/(1+tt.exp(beta[subject_id]*train_x)))
    # you can also just use pm.math.sigmoid() or pm.math.invlogit()
    # which is simpler for a variety of reasons
    #p = pm.math.sigmoid(beta[subject_id] * train_x)
    #p = pm.math.invlogit(beta[subject_id] * train_x)

    # here we feed p to pm.Binomial() to explain train_y_discrete and n_pennies
    obs = pm.Binomial('obs', p=p, n=n_pennies, observed=train_y_discrete)
    trace_new = pm.sample(2000, tune=1000, target_accept=.85)