Bayesian Inference with PyMC3. Compilation error.

Question

The following two codes do a simple bayesian inference in python using PyMC3. While the first code for exponential model compiles and run perfectly fine, the second one for a simple ode model, gives an error. I do not understand why one is working and the other is not. Please help.

Code #1

import numpy as np
import pymc3 as pm

def f(a,b,x,c):
    return a * np.exp(b*x)+c


#Generating Data with error
a, b = 5, 0.2
xdata = np.linspace(0, 10, 21)
ydata = f(a, b, xdata,0.5)
yerror = 5 * np.random.rand(len(xdata))
ydata += np.random.normal(0.0, np.sqrt(yerror))


model = pm.Model()
with model:
    alpha = pm.Uniform('alpha', lower=a/2, upper=2*a)
    beta = pm.Uniform('beta', lower=b/2, upper=2*b)
    mu = f(alpha, beta, xdata,0.5)
    Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=ydata)
    trace = pm.sample(100, tune = 50, nchains = 1)

Code #2

import numpy as np
import pymc3 as pm


def solver(I, a, T, dt):
    """Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
    dt = float(dt)           # avoid integer division
    N = int(round(T/dt))     # no of time intervals
    print N
    T = N*dt                 # adjust T to fit time step dt
    u = np.zeros(N+1)           # array of u[n] values
    t = np.linspace(0, T, N+1)  # time mesh

    u[0] = I                 # assign initial condition
    for n in range(0, N):    # n=0,1,...,N-1
        u[n+1] = (1 - a*dt)*u[n]
    return np.ravel(u)

# Generating data
ydata = solver(1,1.7,10,0.1)
yerror = 5 * np.random.rand(101)
ydata += np.random.normal(0.0, np.sqrt(yerror))

model = pm.Model()
with model:
    alpha = pm.Uniform('alpha', lower = 1.0, upper = 2.5)

    mu = solver(1,alpha,10,0.1)

    Y_obs = pm.Normal('Y_obs', mu=mu, sd=yerror, observed=ydata)

    trace = pm.sample(100, nchains=1)

The error is

Traceback (most recent call last):
  File "1.py", line 27, in <module>
    mu = solver(1,alpha,10,0.1)
  File "1.py", line 16, in solver
    u[n+1] = (1 - a*dt)*u[n]
ValueError: setting an array element with a sequence.

Please help.

Answer 1

The error is in this line:

mu = solver(1,alpha,10,0.1)

You are trying to pass alpha as a value, but alpha is a pymc3 distribution. The function solver only works when you provide a number in the second argument.

The code #1 works because this function

def f(a,b,x,c):
    return a * np.exp(b*x)+c

returns a number.