Bayesian Adaptive Regression Splines and MCMC - Some Questions

Hi all,

I have recently been using this implementation (download available) of Bayesian Adaptive Regression Splines (BARS) in order to attempt to model inventory demand as a non-stationary Poisson process. The implementation is based on this paper.

I have a couple of questions which I am hoping someone here may be able to help me with.

Firstly, does anyone know of a similar implementation in python which I may be able to use?

Secondly, a more general question about MCMC I suppose. I am finding that even with long burn-in and sample runs, getting consistent results from this implementation is proving difficult. The algorithm adds, removes and relocates the spline ‘knots’ in the MCMC run. I am finding that in certain cases, even with a long run, parts of the posterior parameter space are not being explored, i.e. the MCMC appears to be getting stuck in a local minimum.

What is an acceptable approach in cases like this, e.g.

  • Undertake a large number of runs and select those which appear to not get stuck in a local minimum?
  • Find an initial seed value which doesn’t get stuck and use that one run?

In broader terms, how can I ensure that the parameter space has been adequately spanned? Is there any kind of objective measure which can be used?

Many thanks,


You might find the implementation from @AustinRochford helpful

For your second question, both solutions are not recommended. If you are observing problem in the MCMC run, such as sampler being stuck, divergent samples, it is an indication that there are some hidden problem of your model: some area of the parameter space contains high curvature, etc.
I would suggest first try more informative priors to cut out those part of the parameter space, then explore different parameterization of your model.

That implementation looks very useful, thanks, and will certainly be a lot easier to tweak than the C code referenced above.

The problem I face is that I am trying to find a universal method of initialising the MCMC runs for around 1,200 different demand patterns (time series) as I do not have the resource to tweak each one. For instance, I have one simple case where there appear to be local minima at 3 and 5 knots. Sometimes the MCMC will have a modal number of knots of 3, sometimes 5, and sometimes it is bimodal (3 and 5). I can see from the BIC values that usually 3 knots gives a lower minimum BIC - and hence a ‘better model’ by some measure - than 5 knots. But with a run in which the modal number of knots is 5, the lowest BIC is given by a model with 8 knots.

Any help with this would be very welcome.

How is the number of knots being model?

From the implementation paper:

"The algorithm randomly selects as a proposal (in the Metropolis sense) either a new knot (a birth step), removal of a knot (a death step), or a relocation of a knot.

Birth and relocation steps begin by randomly selecting a knot, with equal probabilities, from among all knots in the current knot set. If the proposal involves a birth step, then a location for the potential new knot is randomly selected by first randomly selecting an existing knot and then drawing from a beta distribution centered at that knot. This beta distribution is typically quite tight around its center (with spread controlled by an optional user-defined parameter tau) in order to propose knots that are close to existing knots. The same beta distribution is used to propose a knot relocation (the distribution being centered on the knot to be potentially relocated). As usual, the proposal to add, delete, or relocate is evaluated, and possibly accepted, via the Metropolis-Hastings ratio."

There is pseudocode in the same paper.

So I guess you are modelling as a Discrete variable with something like a Discrete Uniform?

The prior on the number of knots can either be Discrete Uniform in a user-supplied interval, Poisson with a user-supplied lambda, or discrete user-defined (arbitrary).

Hi Matt, same problem with me, that is, can’t find consistent knots location even with long iterations, did you find the solution? Thanks