Thurstone Comparative rankings and 2D shape-parameters

Questions
1

Trying to implement the model found here, which involves deriving a value from observations of comparisons between RV’s.

So, there are some number of things we would like to rank, $X_1, X_2, … X_n$. These are assumed to have a measurable “quality” that is distributed like $X_i ~ N(\mu_i, \sigma_i)$ . A set of judges then sample from this quality, and make a judgement on which quality is higher, i.e., they sample from (assuming no correlation) P(X_i > X_j) = CDF(\mu_i - \mu_j, \sigma_{i,j}). Then, given a whole bunch of judgments (i.e. X_1 > X2, X_5 < X_3, etc.) that form a directed graph, we would like to infer the original parameters mu_1, mu_2, …mu_n.

This boils down to a tt.switch model, but am I right in thinking the shape of the observations would need to be N-by-N? This would mean having some CDF that is shape (n,n) to model the likelihood of observations, right?

EDIT: I’ve gone ahead and attempted something that uses slicing, but I’m not sure how to get Theano to work in this way. Assuming I can “transpose” the PyMC3 RV’s as I can in numpy and theano.tensor.col(), something like this would make sense:

# actual item "quality", and pairwise comparisons
n_itm = 10
a = np.random.rand(n_itm,1)
pwc = np.where(a>a.T, 1,0)

# which comparisons were judged
n_comp = 20
comp_obsv = np.zeros((n_comp,2), dtype=int)
for row in range(n_comp):
    comp_obsv[row, :] = np.random.choice(n_itm, size=2, replace=False)

# dummy observations:
y = pwc[comp_obsv[:,0], comp_obsv[:,1]]

# Thurstone "Law of Comparitive Judgement"
with pm.Model() as model:
    theta = pm.Uniform('true_value', 0,10,shape=n_itm)   # quality priors
    sigma = pm.HalfCauchy('noise', 5, shape=n_itm)  # noise priors 
    
    X = pm.Normal('meas_value', theta, sigma, shape=n_itm)  # model the items
    
    # is each X_i better than each X_j?
    comparison = tt.switch(X>tt.transpose(X), 1, 0)
    
    # pick out the comparisons we actually observed
    judged = comparison[comp_obsv[:,0], comp_obsv[:,1]]
    like = pm.Potential('judgement', judged, observed=y)

But, I’m getting the error:

----------------------------------------------------------------------
IndexError                           Traceback (most recent call last)
<ipython-input-137-9298573097a9> in <module>()
     24 
     25     # pick out the comparisons we actually observed
---> 26     judged = comparison[comp_obsv[:,0], comp_obsv[:,1]]
     27     like = pm.Potential('judgement', judged, observed=y)
     28 

~/anaconda3/lib/python3.6/site-packages/theano/tensor/var.py in __getitem__(self, args)
    538                 return self.take(args[axis], axis)
    539             else:
--> 540                 return theano.tensor.subtensor.advanced_subtensor(self, *args)
    541         else:
    542             if numpy.newaxis in args:

~/anaconda3/lib/python3.6/site-packages/theano/gof/op.py in __call__(self, *inputs, **kwargs)
    613         """
    614         return_list = kwargs.pop('return_list', False)
--> 615         node = self.make_node(*inputs, **kwargs)
    616 
    617         if config.compute_test_value != 'off':

~/anaconda3/lib/python3.6/site-packages/theano/tensor/subtensor.py in make_node(self, x, *index)
   2139 
   2140         index = tuple(map(as_index_variable, index))
-> 2141         bcast = adv_index_broadcastable_pattern(x, index)
   2142         return gof.Apply(self,
   2143                          (x,) + index,

~/anaconda3/lib/python3.6/site-packages/theano/tensor/subtensor.py in adv_index_broadcastable_pattern(a, idx)
   2120     # 2 - True = 1; 2 - False = 2
   2121     fakeshape = [2 - bc for bc in a.broadcastable]
-> 2122     retshape = numpy.empty(fakeshape)[newidx].shape
   2123     return tuple([dim == 1 for dim in retshape])
   2124 

IndexError: too many indices for array

indicating that this kind of slicing (which I thought might work as per the Radon Hierarchical pooling model) isn’t going to work the way I expected. Thoughts?

2

Hi @tbsexton, I am not sure I understand completely what you are trying to do, here are just two quick points:

1, The shape error is resulting from the fact that X is only one dimension. Something like this should work:

    theta = pm.Uniform('true_value', 0,10,shape=(n_itm, 1))   # quality priors
    sigma = pm.HalfCauchy('noise', 5, shape=(n_itm, 1))  # noise priors 
    X = pm.Normal('meas_value', theta, sigma, shape=(n_itm, 1))  # model the items

2, You can not use pm.Potential for observed. You need to write a logp using pm.DensityDist.

3

Thanks! This sent me down the right path, and the model seems to work pretty well.

from itertools import product, combinations, permutations

# actual item "quality", and pairwise comparisons
n_itm = 20
a = np.random.rand(n_itm,1)
a = (a-a.mean())/a.std()
pwc = np.where(a>a.T, 1,0)

# which comparisons were judged
n_comp = 200
combs = np.array(list(permutations(range(n_itm), 2)))
comp_obsv = combs[np.random.choice(combs.shape[0], n_comp, replace=False)]
S = np.zeros_like(pwc)

for obs in comp_obsv:
    S[obs[0],obs[1]] += pwc[obs[0],obs[1]]
    
def norm_cdf(z):
    return 0.5 * (1 + tt.erf(z / np.sqrt(2)))

# Thurstone "Law of Comparitive Judgement"
with pm.Model() as model:
    theta = pm.Normal('true_value', 0, sd=1, shape=(n_itm, 1))   # quality priors
    X = pm.Normal('meas_value', theta, sd=.25, shape=(n_itm, 1))  # model the items
    
    # is each X_i better than each X_j?
    comparison = X-tt.transpose(X)    
    logP = pm.DensityDist('judgement', 
                          lambda Sij: tt.sum(Sij*norm_cdf(comparison)),
                         observed=S)

This samples well and converges nicely with NUTS. The idea is to back out some measure of object “quality” and, therefore, object rankings, when you only know a subset of the true pairwise comparisons between objects.

It’s a way to crowd-source item value functions when the sampling source is a group of people, since people tent to be bad at objective valuations and much better at relative comparisons.

I noticed there’s not too many examples (if any?) on the docs using multi-dimensional observations. I might fix this gist up a bit to illustrate the principle, would this be something of use? Since this works so well, these kinds of observations can be really useful for inferring properties of graphs, even dynamic graphs.

1 Like
4

You are right, we are definitely lacking some examples with multi-dimensional observed. The gist looks interesting, and it reminds me of this example https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/dawid-skene.ipynb

Do you think it might be possible to extend the above example with your model into a coherent doc?