Getting shared parameters with MCMC

Questions
1

Sorry for a newbie question… I have a 4x4 data matrix and would like to obtain parameters alpha (4-vector), beta (4-vector) and c using MCMC such that the posterior distribution in data[i,j] is given by 1/(1 + exp(alpha[i] + beta[j] + c)) (sigmoid). So for 16 data cells I have 9 parameters. I’ve followed the basic pymc3 tutorial and tried the following:

basic_model = pm.Model()
with basic_model:

    # Priors for unknown model parameters
    a = pm.Normal('a', mu=10, sd=10, shape=4)
    b = pm.Normal('b', mu=10, sd=10, shape=4)
    c = pm.Normal('c', mu=0, sd=10)

    alpha = np.tile(a, (4, 1)).T
    beta = np.tile(b, (4, 1))
    offset = np.tile(c, (4, 4))

    # Expected values of outcome 
    # (I'm trying to vectorise the operation here, to get a 4x4)
    mu = pm.math.sigmoid(alpha + beta + offset)  # ERROR!

    # Likelihood (sampling distribution) of observations
    # (data also 4x4 matrix)
    Y_obs = pm.Bernoulli('Y_obs', p=mu, observed=data)

    # draw 500 posterior samples
    trace = pm.sample(50000)

The error is:

AsTensorError: (‘Cannot convert [[Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0\n Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0]\n [Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0\n Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0]\n [Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0\n Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0]\n [Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0\n Elemwise{true_div,no_inplace}.0
Elemwise{true_div,no_inplace}.0]] to TensorType’, <class
‘numpy.ndarray’>)

The error element is probably referring to mu which has too many elements (?). In the tutorial* with the following simple distribution, mu when printed is simply Elemwise{add,no_inplace}.0 but in my mu I have 4x4 elements of Elemwise{true_div,no_inplace}.0. Looks like I’m not allowed to tile the parameters and pass it through. Perhaps it is better to do a matrix factorisation method, if so, how can I go about doing this?

*Tutorial example uses:

alpha = pm.Normal('alpha', mu=0, sd=10) 
beta = pm.Normal('beta', mu=0, sd=10, shape=2) 
sigma = pm.HalfNormal('sigma', sd=1) 
mu = alpha + beta[0] * X1 + beta[1] * X2  # mu = Elemwise{add,no_inplace}.0
2

You can not use numpy function on theano tensors. You can either use the reshaping function in theano.tensor module, or create the shape of RV so that it will broadcast to the shape you want:

basic_model = pm.Model()
with basic_model:
    # Priors for unknown model parameters
    a = pm.Normal('a', mu=10, sd=10, shape=(4, 1))
    b = pm.Normal('b', mu=10, sd=10, shape=(1, 4))
    c = pm.Normal('c', mu=0, sd=10)
    test = a+b <-- this will broadcast to (4 x 4)
    mu = test + c <-- scalar c will add to the (4 x 4) elementwise
1 Like
3

This worked perfectly - thanks so much!