# Order statistics in PyMC3

**URL:** <https://discourse.pymc.io/t/order-statistics-in-pymc3/617>\
**Category:** Questions\
**Created:** [November 30, 2017, 5:02pm UTC](https://discourse.pymc.io/t/order-statistics-in-pymc3/617 "2017-11-30T17:02:34Z")\
**Posts on this page:** 1\
**Showing post:** 3

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [November 30, 2017, 6:39pm UTC](https://discourse.pymc.io/t/order-statistics-in-pymc3/617/3 "2017-11-30T18:39:34Z")

</div>

OK, I cooked up something that kind of work:

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

# data
K = 5 # number of items being ranked
J = 20 # number of raters
yreal = np.argsort(np.random.randn(K, 1), axis=0)
print(yreal)
y = np.argsort(yreal + np.random.randn(K, J), axis=0)
print(y)

class Ordered2D(pm.distributions.transforms.ElemwiseTransform):
    name = "ordered"

    def forward(self, x):
        out = tt.zeros(x.shape)
        out = tt.inc_subtensor(out[0,:], x[0,:])
        out = tt.inc_subtensor(out[1:,:], tt.log(x[1:,:] - x[:-1,:]))
        return out
    
    def forward_val(self, x, point=None):
        x, = pm.distributions.distribution.draw_values([x], point=point)
        return self.forward(x)

    def backward(self, y):
        out = tt.zeros(y.shape)
        out = tt.inc_subtensor(out[0,:], y[0,:])
        out = tt.inc_subtensor(out[1:,:], tt.exp(y[1:,:]))
        return tt.cumsum(out, axis=0)

    def jacobian_det(self, y):
        return tt.sum(y[1:,:], axis=0)

# transformed data{
y_argsort = np.argsort(y, axis=0)

with pm.Model():
    mu_hat = pm.Normal('mu_hat', 0, 1, shape=K-1)
    # set first value to 0 to avoid unidentified model
    mu = tt.concatenate([[0.], mu_hat])
    #sd = pm.HalfCauchy('sigma', 1.)
    latent = pm.Normal('latent',
                       mu=mu[y_argsort],
                       sd=1., # using sd does not work yet
                       transform=Ordered2D(), 
                       shape=(K,J),
                       testval=np.repeat(np.arange(K)[:,None], J, axis=1))
                        # There are some problems with Ordered 
                        # right now, you need to give a testval
    trace = pm.sample()
pm.traceplot(trace, varnames=['mu_hat']);

print('The true ranking is: ')
print(yreal.flatten())
print('The Latent mean is: ')
latentmu = np.hstack(([0], pm.summary(trace, varnames=['mu_hat'])['mean'].values))
print(np.round(latentmu,2))
print('The estimated ranking is: ')
print(np.argsort(latentmu))

```

The Jacobian of the Ordered transformation might not be correct (awaiting for @aseyboldt to confirm as he mentioned to me once).

---

_[View the full topic](https://discourse.pymc.io/t/order-statistics-in-pymc3/617)._