Multi-level multi-variable model with arbitrary and updatable pdf

I am writing a simulated annealing app using pymc3 at each annealing step. Posteriors from each pymc3 run are used as priors for the next, and the “temperature” (assumed noise level) is reduced.

I am using a spline function to calculate my PDF from a uniform distribution over (0,1), and I update the splines after each run. My problem is that I am using several multi-level variables where the splines are distinct for each level, and I cannot figure out how to specify which spline to use for each variable level.

Below is a minimal example of what I am trying to do:

#Define Functions
import numpy as np
import pymc3 as pm
import theano as th
import theano.tensor as tt
import matplotlib.pyplot as plt

def cubic_spline(x, y):
  Creates a spline between points. Interpolations are calculated with Δx from next smaller input point
  x  : list of floats
  y  : list of floats

  list of list of floats
  n = len(x) - 1
  h = [x[i+1]-x[i] for i in range(n)]
  al = [3*((y[i+1]-y[i])/h[i] - (y[i]-y[i-1])/h[i-1]) for i in range(1,n)]
  l = [1] * (n+1)
  u = [0] * (n+1)
  z = [0] * (n+1)
  for i in range(1, n):
    l[i] = 2*(x[i+1]-x[i-1]) - h[i-1]*u[i-1]
    u[i] = h[i]/l[i]
    z[i] = (al[i] - h[i-1]*z[i-1])/l[i]

  b = [0] * (n+1)
  c = [0] * (n+1)
  d = [0] * (n+1)
  for i in range(n-1, -1, -1):    #for i in reversed(range(n)):
    c[i] = z[i] - u[i]*c[i+1]
    b[i] = (y[i+1]-y[i])/h[i] - h[i]*(c[i+1] + 2*c[i])/3
    d[i] = (c[i+1]-c[i])/(3*h[i])
  return [y, b, c, d, x]

def get_pdf(x,a):
    Returns interpolated y value given x input, a is spline created with cubic_spline 
    y = th.shared(np.asarray(a[0]))
    b = th.shared(np.asarray(a[1]))
    c = th.shared(np.asarray(a[2]))
    d = th.shared(np.asarray(a[3]))
    xs = th.shared(np.asarray(a[4]))
    return y[i]+b[i]*(x-xs[i])+c[i]*(x-xs[i])**2+d[i]*(x-xs[i])**3 

def fx(x,q,a,p,ra):
    Returns y values
    return (a+x*p[q]+ra*np.random.random())

def get_pr0(y):
    Returns initial guess PDF for pr[0] values from uniform inputs on (0,1)
    return .01+3.*y

def get_pr1(y):
    Returns initial guess PDF for pr[1] values from uniform inputs on (0,1)
    return .01+10.*y

def get_am(y):
    Returns initial guess PDF for am values from uniform inputs on (0,1)
    return y

#Make Initial Splines



# Generate Dummy Data with constant intercept (am) and slope changing from ps[0] to ps[1] midway through dataset
qs=np.linspace(0,Npts-1,Npts, dtype=int)
for i in range(0,Npts):
    qs[i]=int(2.*i/Npts) # an array to index ps 


with pm.Model() as model:
    am = get_pdf(pm.Beta('am', alpha=1,beta=1),spls[0])
    pr = get_pdf(pm.Beta('pr', alpha=1,beta=1,shape=2),spls[1+qs])    
    eps = Tmpr
    y_pred = pm.Normal('y_pred', mu=mu, sd=eps, observed=ys)
    step = pm.NUTS()
    trace = pm.sample(NSamples,njobs=NJobs,tune=1000)


#Update splines for next iteration, using trace elements (code not shown)

Here I have attempted to index the spline using the indexing array qs, but, of course, I get an error trying to index a list with a numpy array, which doesn’t work, but the problem is that I don’t know what attribute to use to indicate which spline to assign to each parameter level.

This question is similar to Multi-level multi-variables model, in which junpenglao suggested using matrix multiplication shown here, but I’m not sure how I’d implement or update that iteratively for a cubic spline-defined PDF

1 Like

Hmmm, we have a doc which doing a similar thing (smoothing the posterior and use as prior): Have you check it out?

Somehow it seems this notebook is not indexed on the website…