non-linear regression with a parameter vector with pymc

Question

I currently use scipy.optimize.minimize and scipy.optimize.leastsq to perform non-linear regression on my datasets. I would like to use PyMC(3) to investigate the posteriors for all the parameters involved in the fitting procedure. I came across this previous answer on SO.

This is a pretty good example to have available, most of the other examples I saw were for linear regressions. However, the example is not entirely suitable for my purposes. My model has a variable number of parameters, of which I would be fitting a subset. This subset would normally be in the range of 1 to 20 parameters, but sometimes more. With the scipy minimizers those varying parameters are delivered to the cost function in the form of a 1D np.ndarray, p, e.g.

def chi2(p, *args):
    xdata = args[0]
    return p[0] + xdata * p[1] + ........

In the link given above the @pymc.deterministic decorated gauss function has keyword arguments. This is impractical for me, as the same code block needs to deal with varying (and sizeable) numbers of parameters. Is there any way of supplying a vector of parameters instead? I would also have to supply a list of priors for each of the parameters. However, I have a list of lower and upper bounds for each parameter [(min, max)...], so that wouldn't be a problem.

Answer 1

For a set of parameters, it is often best to use a vector and take the dot product:

@deterministic
def theta(beta=beta):
    return beta.dot(X)

or if you want an explicit baseline mean:

@deterministic
def theta(mu=mu, beta=beta):
    return mu + beta.dot(X)

(this is all PyMC 2.3)

Answer 2

Since you are using standard optimisation routine, you may be able to express the function you are minimising as log-likelihood. If it is expressed as log-likelihood and all you want to do is to explore the posterior distributions of your parameters, you may find emcee easier to use. In my case, I was actually minimising a log-likelihood before I started looking into mcmc approach.

from scipy.optimize import minimize

bnds = ((.0001,20),(.0001,20),(.0001,20),(.0001,20))

solution = minimize(my_function_ll, np.array([1.1,1.2,1.3,1.4]), 
                    args=(my_data),
                    bounds=bnds )

All I had to do, is to plug my_function_ll() into emcee like this:

import emcee

# for reproducible results
np.random.seed(0)

ndim = 4  # number of parameters in the model
nwalkers = 10  # number of MCMC walkers
nburn = 1000  # "burn-in" period to let chains stabilize
nsteps = 10000  # number of MCMC steps to take

# we'll start at random locations between 0 and 2
starting_guesses = 2*np.random.rand(nwalkers, ndim)

def log_prior(theta):
    x1,x2,x3,x4 = theta
    # here you can specify boundary. In my case I was using 0 - 20
    if 0.0001 < x1 < 20 and 0.0001 < x2 < 20 and 0.0001 < x3 < 20 and 0.0001 < x4 < 20:
        return 0.0
    return -np.inf

def log_posterior(theta, observation_array):
    lp = log_prior(theta)
    if not np.isfinite(lp):
        return -np.inf    
    return log_prior(theta) + my_function_ll(theta, observation_array) 

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_posterior, 
                                args=[my_data])
sampler.run_mcmc(starting_guesses, nsteps)

sample = sampler.chain[:, nburn:, :].reshape(-1, 4)

Now you can compare MCMC and MLE results:

print "MLE: ", solution.x
print "MCMC: ", sample.mean(0)

Answer 3

The following link contained a lot of information that I was able to use to get going.

https://groups.google.com/forum/#!msg/pymc/WUg0tZjNAL8/jsVRFe1DMT4J