Pymc3- inside pm.Model() trying to pass dictionary, NOT working

Questions
1

The following code is a adaptation of a code from web(I don’t remember the source). It does a simple Bayesian Inference using PyMC3 for a Lotka-Volterra model.

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import pymc3 as pm
import theano.tensor as T
import os, sys
from theano.compile.ops import as_op

#Lotka volterra solution for given input parameters
def lv(a, b, c, I):
    # growth rates
    def dX_dt(X, t=0):
        #RHS of the Lotka Volterra equations
        return np.array([ a*X[0] -   b*X[0]*X[1] , 
                        -c*X[1] + I*X[0]*X[1] ])
    # Initial conditions
    X0 = np.array([2, 3])
    # Solving the model for given parameter set
    X, infodict = odeint(dX_dt, X0, times, full_output=True)
    rabbits, foxes = X.T
    
    return rabbits, foxes

times = np.linspace(0, 150,  150)

#Parameter values for generating synthetic data
a0=2
b0=8
c0=2
I0=4

#Generating Synthetic Data with fixed parameters
rdata, fdata = lv(a0, b0, c0, I0)

#Adding some errors to the generated data
rerror = np.random.rand(len(times))
ferror = np.random.rand(len(times))
rdata += 0.5*np.random.normal(0.0, rerror**2)
fdata += 0.5*np.random.normal(0.0, ferror**2)

#Final Data - ydata is the observed values and yerror are the corresponding errorbars.
ydata = np.ravel((rdata, fdata))
yerror = np.ravel((rerror, ferror))

#This part I don't understand clearly
@as_op(itypes=[T.dscalar, T.dscalar, T.dscalar, T.dscalar], otypes=[T.dvector])

def solution(a, b, c, I):
    v, w = lv(a=a, b=b, c=c, I=I)
    return np.ravel((v, w))
    
#MCMC Sampling

nsample = 5000
burn = 2000

model = pm.Model()
with model:
	#Priors
	a = pm.Uniform('a', lower=a0/2, upper=a0+a0/2)
	b = pm.Uniform('b', lower=b0/2, upper=b0+b0/2)
	c = pm.Uniform('c', lower=c0/2, upper=c0+c0/2)
	I = pm.Uniform('I', lower=I0/2, upper=I0+I0/2)
	#Solving the system
	mu = solution(a, b, c, I)
	#Likelihood
	Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=ydata)
	trace = pm.sample(nsample, pm.Metropolis(), tune = burn, nchains = 1)

First I generate the predator and prey population data for given period of time by solving the differential equations using known values of parameters. Then I add some errors and errorbars to this to make it look realistic.
I use PyMC3 to draw th posterior of the parameters. This works fine.
But now I am using a much more complex model and has a lot of parameters. So I want to do it as compactly as possible. I would like to pass a dictionary of parameters to the solution function and use it to solve the model. I can pass the parameters one by one but it is a lot of hard coding. Is there any workaround?

2

Hi!

Why must your parameters be in a dictionary? For instance, the following code will work:

with pm.Model() as model:
    a = pm.Uniform('a', lower=0, upper=1)
    b = pm.Uniform('a', lower=0, upper=1)
    sol = pm.Deterministic('sol', a+b)  # Or whatever sol you need
    err = pm.HalfNormal('err', sd=1)  # Or whatever err you need
    likelihood = pm.Normal('likelihood', mu=sol, sd=err, observed=data)

    pm.sample()

Then, you can assign parameters like a and b priors, line by line. Afterwards, you can access your parameters like so: trace['a']

Perhaps I’m misunderstanding how you want to use PyMC3?

3

That is not what I intended. Sorry if the question was obscure. I have modified it. Please see it and help.

4

Hello!

Sorry for misunderstanding the first time around. I’m still a bit confused by what you mean when you say you’re “using a much more complex model and has a lot of parameters”. Either you mean

a) your model has a lot of the similar parameters (i.e. there may be five different a parameters, which will be assigned the same prior but may end up having different posteriors), or
b) your model has more, new parameters (i.e. there are completely different parameters which will be assigned different priors)

If a), then you should use the shape keyword. e.g. x = pm.Uniform('x', lower=0, upper=1, shape=[5, 5]) will give you a 5x5 matrix of uniform random variables.
If b), then there isn’t any good way to save on typing. Nor should there be: if your priors are different, they must be explicitly specified.

Hope this helps!