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
Parameters
----------
x : list of floats
y : list of floats
Returns
-------
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)]
al.insert(0,0)
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
"""
i=tt.cast(np.floor((x-a[4][0])/(a[4][1]-a[4][0])),'int64')
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
xi=np.linspace(0.,1.,100)
yi=get_am(xi)
sp0=cubic_spline(xi,yi)
sp1=cubic_spline(xi,get_pr0(xi))
sp2=cubic_spline(xi,get_pr1(xi))
spls=[sp0,sp1,sp2]
Tmpr=1
NSamples=1000
NJobs=1
Npts=10000
# Generate Dummy Data with constant intercept (am) and slope changing from ps[0] to ps[1] midway through dataset
xs=np.linspace(0,20,Npts)
ps=np.linspace(1.,2.,2,dtype=float)
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
ys=np.zeros((Npts))
ys=fx(xs,qs,.3,ps,0.)
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
mu=fx(xs,qs,am,pr,0.)
y_pred = pm.Normal('y_pred', mu=mu, sd=eps, observed=ys)
step = pm.NUTS()
trace = pm.sample(NSamples,njobs=NJobs,tune=1000)
pm.traceplot(trace)
pm.autocorrplot(trace);
plt.show()
#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