CODEWITHSUNDEEP
pythonpymc3
Olivier Ma
Olivier Ma

Reputation: 1309

Pymc3: how to defined an ordered vector of parameters

I'm trying to do an ordered logistic regression, one parameter of which is an ordered vector of cutpoints. I'm not sure how I can define them.

One very stupid way I have come up with is to just define each component of the vector manually, using one as the other's bound:

with pm.Model() as bound_model:
    a = pm.Normal('a', mu=0, sd=10)
    BoundedB = pm.Bound(pm.Normal, upper=a)
    b = BoundedB('b', mu=0, sd=10)
    BoundedC = pm.Bound(pm.Normal, upper=b)
    c = BoundedC('c', mu=0, sd=10)

    bound_trace = pm.sample(1000)

This is hardly efficient and I'm not sure whether they will work as expected. Is there better way to do this?

Upvotes: 0

Views: 848

Answers (1)

aseyboldt
aseyboldt

Reputation: 1090

This is a missing feature in pymc3 I guess. I might write a pull request, but in the meantime you can use something like this:

class Ordered(pymc3.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 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)

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

And use it like this:

N = 10
with pm.Model() as model:
    pm.Normal('y', mu=np.zeros(N), sd=1, shape=N, 
              transform=Ordered(), testval=np.arange(N))

Edit: A very short explanation about what is happening here:

We define a map from $R^n$ to the set of ordered sequences by

f(x_1) = x_1,\quad f(x_i) = f(x_{i - 1}) + exp(x_i)

Since this is a nice bijective function, we can compute the probability density on $R^n$ by

P_{R^n}(x) = P_{ordered}(f(x)) \cdot |J_{f(x)}|

where J is the jacobian of the transformation.

The sampler will only see the unconstrained values. This is pretty much how Bound is implemented in the first place.

If you want more details, you can have a look at the stan manual. It contains a nice description of these transforms, and the math is the same for pymc3 and stan.

Upvotes: 4

Related Questions