Setting subtensors within the model code block

Questions
1

Hello, this is my first post on discourse - thank you for the knowledge and assistance!

I’m an environmental engineer working with isotope tracers in sediment (lead210 and cesium137). There are some other Bayesian approaches for isotope dating but not exactly what I’m looking for. My desired model is to have a concentration (distribution) for each year period; a mixing depth (based on the same annual layers, e.g., the upper 10 years of sediment mix); and radioactive degradation over time. I can generate data in numpy/ pandas and with theano tensors, but I fall on my face within the pymc3 model.

To start, I’m trying to get (just) one year of degration/ deposition/ and mixing to work in the model, but I’m not able to get it to work. Here is the code (hopefully it formats):

import pymc3 as pm
import numpy as np
import theano.tensor as tt

model = pm.Model()

#priors
w = 0.1                             # g/cm2/yr
Phi = 10                           # Bq/cm2/yr
simulation_time = 100     # years
layers = 50                      # layers in the sediment core
hl = 22.3                         # half life
lam  = np.log(2)/hl          # decay constant (yr)
an_deg = np.exp(-lam)   # annual degradation

with model:
    # Prior for Phi
    Phi = pm.Gamma('Phi', mu = Phi, sigma = Phi/3, shape = simulation_time )      # Bq/cm2/yr
    
    #populate with the first value all the way down [could degrade .   ]
    At = pm.Deterministic('At', tt.zeros((layers)).fill(Phi[0]/w))
    # annual degradation
    At = At*an_deg
    #shift down
    At = tt.set_subtensor(At[1:], At[0:-1])
    # deposit
    At = tt.set_subtensor(At[0], Phi[1]/w)  # just trying to put any value in the array
    # mix
    At = tt.set_subtensor(At[0:ma], At[0:ma].mean())

sams  = 50
pp = pm.sample_prior_predictive( samples=sams, model=model)

pp['At'][1]

The variable At is not being modified by set_subtensor and remains uniform for any sample draw. I would like to know how to update it.

From there any insight on how to loop through 100 such iterations would be welcome. I’m struggling to understand how this notebook applies.

Let me know if you need more information, thanks!

2

It looks to me like your problem is in the line At = pm.Deterministic('At', tt.zeros((layers)).fill(Phi[0]/w)) - the values are all uniform from there, I believe because you’re using Phi[0]

I got a non-uniform set of samples using

with pm.Model() as model:
    Phi_ = pm.Gamma('Phi', mu = Phi, sigma = Phi/3, shape=(simulation_time,))
    At = pm.Deterministic('At', tt.zeros(simulation_time).fill(Phi_/w))
    pp = pm.sample_prior_predictive(50)

however I had to replace your layers value with simulation time to make the dimensions match

Minor convenience note - ma in your example isn’t defined, and Phi is defined twice…

3

Thanks for looking into this! My intent is to start At uniform and modify during the algorithm. I found that I could get it to work by manipulating At prior to making it a Deterministic variable. Seems strange, but it worked as intended.

I’m currently struggling with scan but need to bang my head on the wall a bit more before reaching out for help. I’ll make sure the variables are clearer in any future post.

Thanks again.