Question: Specifying two deterministic self-referential timeseries

Questions
1

Hello all! I’m a PyMC novice, so sorry if this is the wrong place to ask this. Please feel free to direct me to the correct place.

I’m looking for some guidance in setting up a model-fitting procedure I’m working on. The model basically have 5 stochastic priors (U, R, D, F, f, and A), and two deterministic “time-series” variables (ui and Rn). ui and Rn depend on their initial conditions (U and R respectively). In addition, Rn[idx] depends on ui[idx]. I’ve been trying to make this work for a few days now, and I’m at a place where I need to ask for help.

I’ve tried a whole host of things that haven’t worked (i.e. days of error messages). I’ll spare you all that, and just show you what I have right now. I’m sure there’s a sensible way to accomplish this, and I’m hoping if I just share what I have, someone will be able to tell me a better way to handle it:

The function to calculate the updates to the ui and Rn variables is:

def calc_both(ui, R, U, f, dT, F, D):

ui_main_term = ui + f * (1 - ui) - U
ui_exponential = np.exp(-dT/F)
ui_full_term = U + (ui_main_term * ui_exponential)

Rn_main_term = 1 - (R * (1 - ui))
Rn_exponential = np.exp(-(dT/D))
Rn_full_term = 1 - (Rn_main_term * Rn_exponential)
    
return ui_full_term, Rn_full_term

And the model itself is specified as:

with pm.Model() as TM_model:

# Set constants
dT = 0.1
stim_n = len(data)
R = 1.0

# Set priors
U = pm.Uniform('U', lower=0, upper=1)
f = pm.Uniform('f', lower=0, upper=1)
D = pm.Uniform('D', lower=0, upper=2)
F = pm.Uniform('F', lower=0, upper=2)
A = pm.Uniform('A', lower=20, upper=1000000)
ui = U

# Calculate dynamics
result, updates = tn.scan(fn=calc_both, outputs_info=[ui, R], non_sequences=(U, f, dT, F, D), n_steps=stim_n)
power = tn.function(inputs=[ui, R, U, f, dT, F, D], outputs=result, updates=updates, on_unused_input='warn')
ui_Rn = power(ui, R, U, f, dT, F, D)
ui = ui_Rn[0]
Rn = ui_Rn[1]

# Predict PSCs
p_psc = A * ui * Rn

# Define Likelihood
p_data = pm.Normal('p_data', mu=p_psc, observed=data)

This particular example fails when setting up the Theano Scan function. I’ve tried implementing the loop in pure python, but that fails as well. Any advice would be appreciated!

TL;DR: I’m looking for a working way to specify two deterministic variables (ui, and Rn). They both depend on their initial condition (U and R respectively) and Rn depends on ui. Any help, guidance, or redirection is appreciated!

2

I think you are on the right path - you just need to make sure the tt.scan part is implemented correctly (which is the most challenging part!). Have a look at eg: https://docs.pymc.io/notebooks/PyMC3_tips_and_heuristic.html#PyMC3-implementation-using-theano.scan about how to implement a tt.scan and compare with python for loop

3

Thanks for the quick feedback! I’d certainly say it’s challenging, particularly since i’m new to Theano in general.

So I’ve made a small amount of progress and the theano.scan it works when I run the test outside the PyMC3 model. To accomplish this, instead of explicitly specifying my constants I replaced them with tt.fscalar() values. But when I test it within a pm.Model() it fails with error:

AttributeError: stim_n has no test value

It seems PyMC always sets theano.config.compute_test_value to ‘raise’ when i instantiate the model. I’ve tried running “tn.config.compute_test_value = ‘off’” at the start of the model, but this causes the error:

AttributeError: ‘scratchpad’ object has no attribute ‘test_value’

So I’m guessing that’s a bad idea? Through searching I’ve found several variants of this issue, and they all seem to happen sometimes to some people and not others. Maybe due to broken packages? As per this issue, I just ran theano.test() and had 261 errors :sweat_smile:. I guess it’s time to try a fresh install!

I’ll update with my progress here in case someone finds it helpful in the future. Again, if you know of anything obvious I’m missing please let me know!

4

If you are getting:

One of the simple solution is to do something like mytensor.tag.test_value = 1.