Pymc3: как определить упорядоченный вектор параметров

Question

Я пытаюсь сделать упорядоченную логистическую регрессию, одним из параметров которой является упорядоченный вектор точек среза. Я не уверен, как я могу их определить.

Один очень глупый способ, который я придумал, - просто определить каждый компонент вектора вручную, используя один в качестве границы другого:

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)

Это вряд ли эффективно, и я не уверен, будут ли они работать так, как ожидалось. Есть ли лучший способ сделать это?

Answer 1

Это отсутствующая функция в pymc3. Я мог бы написать запрос на получение, но в то же время вы можете использовать что-то вроде этого:

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:])

И используйте это так:

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

Изменить: очень краткое объяснение того, что здесь происходит:

Определим отображение из $R^n$ на множество упорядоченных последовательностей

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

Поскольку это хорошая биективная функция, мы можем вычислить плотность вероятности на $R^n$

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

где J - якобиан преобразования.

Сэмплер будет видеть только неограниченные значения. Это довольно много, как Bound реализуется в первую очередь.

Если вы хотите больше подробностей, вы можете взглянуть на руководство по стану. Он содержит хорошее описание этих преобразований, и математика одинакова для pymc3 и stan.