Search code examples
pythonbayesianpymc

PyMC, deterministic nodes in loops


I'm a bit new to Python and PyMC, and making rapid progress. But I'm just confused about the use of setting deterministic values of a 2D matrix. I have a model below, that I cannot get to parse correctly. The problem relates to setting the value theta in the model.

import numpy as np
import pymc

define known variables

N = 2
T = 10
tau = 1

define model... which I cannot get to parse correctly. It's the allocation of theta that I'm having trouble with. The aim to to get samples of D and x. Theta is just an intermediate variable, but I need to keep it as it's used in more complex variations of the model.

def NAFCgenerator():

    D = np.empty(T, dtype=object)
    theta = np.empty([N,T], dtype=object)
    x = np.empty([N,T], dtype=object)

    # true location of signal
    for t in range(T):
        D[t] = pymc.DiscreteUniform('D_%i' % t, lower=0, upper=N-1)

    for t in range(T):
        for n in range(N):
            @pymc.deterministic(plot=False)
            def temp_theta(dt=D[t], n=n):
                return dt==n
            theta[n,t] = temp_theta

            x[n,t] = pymc.Normal('x_%i,%i' % (n,t),
                               mu=theta[n,t], tau=tau)

    return locals()

** EDIT **

Explicit indexing is useful for me as I'm learning both PyMC and Python. But it seems that extracting MCMC samples is a bit clunky, e.g.

D0values = pymc_generator.trace('D_0')[:]

But I am probably missing something. But did I managed to get a vectorised version working

# Approach 1b - actually quite promising
def NAFCgenerator():

# NOTE TO SELF. It's important to declare these as objects
D = np.empty(T, dtype=object)
theta = np.empty([N,T], dtype=object)
x = np.empty([N,T], dtype=object)

# true location of signal
D = pymc.Categorical('D', spatial_prior, size=T)

# displayed stimuli
@pymc.deterministic(plot=False)
def theta(D=D):
    theta = np.zeros([N,T])
    theta[0,D==0]=1
    theta[1,D==1]=1
    return theta

#for n in range(N):    
x = pymc.Normal('x', mu=theta, tau=tau)

return locals()

Which seems easier to get at MCMC samples using this for example

Dvalues = pymc_generator.trace('D')[:]

Solution

  • In PyMC2, when creating deterministic nodes with decorators, the default is to take the node name from the function name. The solution is simple: specify the node name as a parameter for the decorator.

            @pymc.deterministic(name='temp_theta_%d_%d'%(t,n), plot=False)
            def temp_theta(dt=D[t], n=n):
                return dt==n
            theta[n,t] = temp_theta
    

    Here is a notebook that puts this in context.