Детерминированные переменные и функция Fortran Scipy в PyMC 3

Я пытаюсь построить простую модель PyMC 3, в которой я оцениваю две точки отсечения и параметр корреляции в скрытой двумерной гауссовой плотности, производя четыре предсказанные вероятности для вектора (полиномиального) подсчета. (Я надеюсь, что в конечном итоге это станет частью более крупной модели, в которой эти и другие параметры оцениваются для ряда скрытых многомерных гауссовских плотностей.)

Итак, я хочу смоделировать точки отсечения cx и cy как обычные случайные величины, а параметр корреляции rho — как масштабированную бета-случайную величину (в качестве примечания, я хотел бы услышать лучший способ справиться с rho — делает ли PyMC 3 усечены нормальные случайные величины, например?). И я хочу использовать функцию mvnun для вычисления предсказанных вероятностей для заданных значений cx, cy и rho. Функция mvnun является частью scipy.stats.mvn, фрагмента скомпилированного кода на Фортране с двумя функциями для вычисления очень точных многомерных нормальных значений CDF.

Если я попытаюсь вставить rho в корреляционную матрицу S или поместить cx или cy в массивы, указывающие пределы интегрирования, я получу следующее:

ValueError: setting an array element with a sequence.

Если я использую фиксированные числовые значения для cx, cy и/или rho, mvnun работает нормально (в блоке 'with model:' или вне его). Я копался, пытаясь понять, почему RV PyMC вызывают эту ошибку, но я в тупике. Я понимаю, что cx, cy и rho являются теано-тензорными переменными, но я не могу понять, что, если вообще что-то, в теано-тензорных переменных может вызвать эти проблемы.

Есть ли фундаментальная проблема с попыткой вызвать функцию Fortran с PyMC RV в качестве аргументов? Или мой код как-то неисправен? Оба? Что-то совсем другое?

Я новичок в PyMC, и я установил PyMC 3, думая, что это текущая версия (клянусь, примечания о том, что это альфа-версия, не было, когда я устанавливал ее несколько недель назад). Возможно, мне следует установить 2.3 и выяснить, как соединить это с этим?

В любом случае, любые советы о том, как исправить ситуацию, будут высоко оценены.

Вот мой код:

import pymc as pm
import numpy as np
from scipy.stats.mvn import mvnun as mvncdf

counts = np.array([100,35,45,20])
N = np.sum(counts)

def scaleBeta(x):
    return x*1.98 - 0.99

model = pm.Model()

with model:

    cx = pm.Normal('Cx', mu=0., tau=1.)
    cy = pm.Normal('Cy', mu=0., tau=1.)

    mu = np.array([0., 0.])

    rho_beta = pm.Beta('rho_beta', alpha=2, beta=2)
    rho = pm.Deterministic('rho',scaleBeta(rho_beta))
    S = np.array([1, rho, rho, 1]).reshape(2,2)

    low_aa = np.array([-50.,-50.])
    upp_aa = np.array([cx, cy])
    low_ab = np.array([-50., cy])
    upp_ab = np.array([cx, 50.])
    low_ba = np.array([cx, -50.])
    upp_ba = np.array([50., cy])
    low_bb = np.array([cx, cy])
    upp_bb = np.array([50., 50.])

    p_aa,i = mvncdf(lower=low_aa,upper=upp_aa,means=mu,covar=S)
    p_ab,i = mvncdf(lower=low_ab,upper=upp_ab,means=mu,covar=S)
    p_ba,i = mvncdf(lower=low_ba,upper=upp_ba,means=mu,covar=S)
    p_bb,i = mvncdf(lower=low_bb,upper=upp_bb,means=mu,covar=S)

    pv = np.array([p_aa,p_ab,p_ba,p_bb])

    cv = pm.Multinomial('counts',n=N,p=pv,observed=counts)