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.

The number of parameters should not change from iteration to iteration, should it? Can you elaborate on this? — Chris Fonnesbeck, Nov 30 '14 at 23:45
Chris, I've trying to bolt on pymc to the lmfit-py project. The idea was to use pymc to perform a Bayesian analysis for parameters in a curvefitting scenario. In lmfit one creates an objective function that returns the residuals for the curvefit. I've since worked out how to use pymc in this case, it's great. Essentially I create a `@pm.observed likelihood` function, that returns a `pm.normal_like` distribution based on the users objective function. In lmfit your curvefit may have M parameters, of which you are allowing N of them to vary, 0 < N <= M. — Andrew Nelson, Dec 30 '14 at 23:57

score 3 · Answer 1 · answered Dec 09 '14 at 23:13

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)

score 2 · Answer 2 · answered Nov 30 '14 at 23:47

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)

score 0 · Answer 3 · answered Dec 31 '14 at 00:00

0

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

answered Dec 31 '14 at 00:00

Andrew Nelson

460
3
11

non-linear regression with a parameter vector with pymc

3 Answers3