How to build the MMSB model with pymc3?

Questions
1

Here is the model I want to describe in pymc3.

Untitled
Untitled.jpg841×351 35.6 KB

I wrote the following codes, but they do not work.

with pm.Model() as model:   
    bernoulli_parameters = np.empty([N, N], dtype=object)
    for n in range(N):  
        for m in range(N):
            if (z[n,m]==z[n,m]):
                bernoulli_parameters[n,m] = myBeta[z[n,m]]
            else:
                bernoulli_parameters[n,m] = epsilon
    bernoulli_parameters_vector=bernoulli_parameters.reshape(1,N*N)    
    D = pm.Bernoulli('D',                                 
            p=bernoulli_parameters_vector, observed=Topology_Data_Vector)

I think it is due to the following 2 reasons：

  1. As both z[n,m] and z[m,n] are hidden variables, I found I can not compare z[n,m] with z[m,n]. (z[n,m] == z[m,n] is always false).
  2. I can not use hidden variable z[m,n] as the index of myBeta, namely myBeta[z[m,n]] also does not work.
2

If myBeta and epsilon is not a tensor, you can move them outside of the with pm.Model()... - it’s much easier to check the correctness of the output and much faster.

In general, you can get rid of the two forloop with some matrix tricks (usually that is way faster than for loop). For example, in your case you can try something like:

import theano.tensor as tt
z = np.random.randint(0, 5, size=(5, 5))
eps = .1
with pm.Model() as m:
    myBeta = pm.Beta('b', 1., 1., shape=5)
    pD = tt.switch(tt.eq(z, z.T), myBeta[z], eps)
    y = pm.Bernoulli('D', p=pD, observed=np.random.rand(*z.shape)>.3)

The nice thing of the code above is that, it works with z being either a numpy array or a theano tensor.

1 Like
3

Thanks so much for your reply.

In the model, z follows a categorical distribution. But the following code still does not work.
(For your convenience, you can download the codes here:etection.py (2.0 KB))

with pm.Model() as model:   
 
    myBeta = np.eye(N)*0.8
    myBeta = myBeta + np.ones([N,N])*0.2-np.eye(N)*0.2   
    
    pi = pm.Dirichlet('pi', a=alpha*np.ones(K),shape=[N,K])    
   
    z = np.empty([N,N], dtype=object)    
    for node1 in range(N):        
            z[node1] = pm.Categorical('z_%d' % node1, p=pi[node1],shape=N)    
        
    pD = tt.switch(tt.eq(z, z.T), myBeta[z], eps)
    y = pm.Bernoulli('D', p=pD, observed=Topology_Data)

The error messages are as follows:

> c:\users\twp\downloads\ection.py(64)<module>()    
3--> 64     pD = tt.switch(tt.eq(z, z.T), 
ipdb> myBeta[z], eps)
     65     y = pm.Bernoulli('D', p=pD, observed=Topology_Data)

theano.tensor.var.AsTensorError: ('Cannot convert [[z_0 z_0 z_0 z_0 z_0]\n [z_1 z_1 z_1 z_1 z_1]\n [z_2 z_2 z_2 z_2 z_2]\n [z_3 z_3 z_3 z_3 z_3]\n [z_4 z_4 z_4 z_4 z_4]] to TensorType', <class 'numpy.ndarray'>)
4

In general, PyMC3 operates on theano tensors, so in your code above, you cannot treat z as numpy array.

Vectorizing z into a tensor should do:

with pm.Model() as model:
    myBeta = pm.Beta('b', g_0, h_0, shape=K)
    pi = pm.Dirichlet('pi', a=alpha*np.ones(K),shape=[N,K])
    z = pm.Categorical('z', p=pi, shape=(N, N))
    pD = tt.switch(tt.eq(z, z.T), myBeta[z], eps)
    y = pm.Bernoulli('D', p=pD, observed=Topology_Data)

Note that:
1, You specify the shape of z, so you dont need a for loop to create it.
2, myBeta[z] broadcast shape = (K,) RV myBeta into the same shape as z with shape = (N, N)
3, I did not reshape the Topology_Data, it still has a shape of (N, N)

1 Like