Using stan as an MCMC kernel within a bigger algorithm

[david...@gmail.com]

Thanks to everyone involved for a very polished piece of software!  Got up and running very quickly (and without a background in bugs!).

I have a question about the viability of using Stan as an implementation of HMC/NUTS for an MCMC kernel - but used in combination with other MCMC kernels that I would like to write myself.

There are a few interesting applications such as:

* a way to handle discrete parameters (i.e. do that bit outside Stan) - although sounds like this will be solved shortly anyway...

* infinite dimensional models.

* applying kernels to annealed distributions e.g. annealed importance sampling

* explicitly using tricks like this: http://andrewgelman.com/2011/10/19/an-interweaving-transformation-strategy-for-boosting-mcmc-efficiency/

My guess is that while this is non-trivial - it can be done by "simply" writing R or python wrappers for a bunch of underlying functions.  Any thoughts?  Is this a common request (sorry if I missed it).

Cheers,
David

[david...@gmail.com]

I posted this too quickly... now I see that there are functions like "sampling" that should expose the required parts of Stan....

Still, I am curious to know if there is much demand for this type of approach...

[Ben Goodrich]

On Monday, July 21, 2014 9:40:47 PM UTC-4, david...@gmail.com wrote:

I have a question about the viability of using Stan as an implementation of HMC/NUTS for an MCMC kernel - but used in combination with other MCMC kernels that I would like to write myself.

Any kind of an external subroutine would be challenging. For continuous unknowns, at a minimum, you have to get the gradient information from the external subroutine back into the Stan internals, and there is no way to do that except to write custom C++ code.
 

There are a few interesting applications such as:
* a way to handle discrete parameters (i.e. do that bit outside Stan) - although sounds like this will be solved shortly anyway...

Won't be solved shortly, although it is probably the most common request.
 

* infinite dimensional models.

There was some thread about Dirichlet Processes recently
 

* applying kernels to annealed distributions e.g. annealed importance sampling

Pretty sure this will never be in Stan, but thermodynamic sampling will probably be in Stan 3.x

http://andrewgelman.com/2014/05/20/thermodynamic-monte-carlo/
 

* explicitly using tricks like this: http://andrewgelman.com/2011/10/19/an-interweaving-transformation-strategy-for-boosting-mcmc-efficiency/

I don't know about that particular trick but the centered vs. non-centered parameterizations is something people have been thinking about

http://arxiv.org/abs/1312.0906

Ben

[david...@gmail.com]

On Tuesday, July 22, 2014 1:50:43 PM UTC+10, Ben Goodrich wrote:

On Monday, July 21, 2014 9:40:47 PM UTC-4, david...@gmail.com wrote:

I have a question about the viability of using Stan as an implementation of HMC/NUTS for an MCMC kernel - but used in combination with other MCMC kernels that I would like to write myself.

Any kind of an external subroutine would be challenging. For continuous unknowns, at a minimum, you have to get the gradient information from the external subroutine back into the Stan internals, and there is no way to do that except to write custom C++ code.

Thanks for the response!

So if I run sampling() and then extract the last element of the chain and pass it back to another call to the init argument in sampling() that won't work? I had almost convinced myself it was that easy....

Cheers,
David

[david...@gmail.com]

To clarify here is some code:

library(rstan)

m <- stan_model(model_code = 'parameters {real y;} model {y ~ normal(0,1);}')

f <- sampling(m, iter = 1,chains=1)
for (s in 1:1000) {
    f <- sampling(m, iter = 1,chains=1,init=list(list('y'=extract(f)$y[1],'lp__'=extract(f)$lp__[1])))
    // apply another kernel to f here (or more likely p(z,f) to update z and pass z back as data.)
    // recompute lp__
}

Is this not valid? 
Thanks,
David

[Ben Goodrich]

On Tuesday, July 22, 2014 12:13:53 AM UTC-4, david...@gmail.com wrote:

To clarify here is some code:

library(rstan)

m <- stan_model(model_code = 'parameters {real y;} model {y ~ normal(0,1);}')

f <- sampling(m, iter = 1,chains=1)
for (s in 1:1000) {
    f <- sampling(m, iter = 1,chains=1,init=list(list('y'=extract(f)$y[1],'lp__'=extract(f)$lp__[1])))
    // apply another kernel to f here (or more likely p(z,f) to update z and pass z back as data.)
    // recompute lp__
}

Is this not valid?

If the recomputed lp__ depends on the continuous unknown y, then HMC would need to know the derivative of that with respect to y. And that would basically entail working out the derivative analytically and writing custom C++. Moreover, Stan currently has a really hard distinction between known data (that never changes) and unknown parameters (that it updates). So, there is not a mechanism by which you can pass "data" that is really an unknown that gets updated some other way.

Ben

[Michael Betancourt]

Aside from the technical difficulties that Ben notes, it’s probably not going to be a good

idea in practice.  For example, see the references he noted for why simulated annealing

(and simulated importance sampling) do not scale well and why interweaving parameterizations

won’t help much for a proper joint sampler (what good taste in papers, Ben!).

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[david...@gmail.com]

Thanks Ben,

I wrote a slightly more elaborate script which seems to confirm your point...  it seems something is being forgotten or reset in the Markov chain and it doesn't have detailed balance.. (probably as you say the derivatives).  It still puzzles me why this is necessary... (as Stan can compute these - and they cannot be passed to sampling()).

Here I hoped to get around the limitations by making the Stan kernel apply to p(y|mu)  and then sampling mu straight from p(mu) - but it doesn't seem to work (although poor mixing might also be the problem).

I hear what you say about the hard distinction, between data and parameters - but what that means in practice still puzzles me....   The way I see it if everything is collected as needed from sample() and the passed back to it - I don't see why this can't work.  Maybe momentum / derivatives are the key problem... (although I thought momentum was an independent draw)

m <- stan_model(model_code = 'data {real mu;} parameters {real y;} model {y ~ normal(mu,1);}')

innersamp <- 1
numsamples <- 10000
f <- sampling(m, ,data=list('mu'=mu),iter = innersamp,warmup=0,chains=1)
samp <- matrix(0,1,numsamples)

for (ii in 1:numsamples) {
    loglike <- -0.5*(extract(f)$y-mu)^2
    f <- sampling(m, iter = innersamp,warmup=0,chains=1,data=list('mu'=mu),init=list(list('y'=extract(f)$y[innersamp],'lp__'=loglike)))
    mu <- rnorm(1,0,1)
    samp[ii] <- extract(f)$y[innersamp];
}

print(apply(samp,1,sd))# should be root 2  (seems to be lower)
# samp[1:numsamples-1]==samp[2:numsamples] # check accept/reject

[david...@gmail.com]

On Tuesday, July 22, 2014 6:22:32 PM UTC+10, Michael Betancourt wrote:

Aside from the technical difficulties that Ben notes, it’s probably not going to be a good

idea in practice.  For example, see the references he noted for why simulated annealing

(and simulated importance sampling) do not scale well and why interweaving parameterizations

won’t help much for a proper joint sampler (what good taste in papers, Ben!).

My hope was primarily to deal with models changing dimension and with discrete latent variables...   but I put in a couple of other suggestions ...  (I agree the annealing bit is particularly questionable)...  On the other hand Mark Girolami actually recommends interweaving here: http://www.dcs.gla.ac.uk/publications/PAPERS/9382/tech_rep_lgms.pdf .  Iain Murray argues running HMC on the full posterior for latent variable Gaussian process hierarchical  models doesn't work so well - but I must admit I don't understand why.

It is fairly common that a published algorithm uses HMC as a kernel in combination with standard Metropolis-Hastings or Gibbs kernels...  it would be surprising if this were _never_ a good idea....   It would really be wonderful to be able to use Stan to make the HMC part easier...  but it looks like that is at least not easy right now...  although I am still a little puzzled as to why...

Thanks again, for the feed back and for Stan!

 

[Michael Betancourt]

My hope was primarily to deal with models changing dimension and with discrete latent variables...   but I put in a couple of other suggestions ...  (I agree the annealing bit is particularly questionable)...  On the other hand Mark Girolami actually recommends interweaving here: http://www.dcs.gla.ac.uk/publications/PAPERS/9382/tech_rep_lgms.pdf .  Iain Murray argues running HMC on the full posterior for latent variable Gaussian process hierarchical  models doesn't work so well - but I must admit I don't understand why.

It is fairly common that a published algorithm uses HMC as a kernel in combination with standard Metropolis-Hastings or Gibbs kernels...  it would be surprising if this were _never_ a good idea....   It would really be wonderful to be able to use Stan to make the HMC part easier...  but it looks like that is at least not easy right now...  although I am still a little puzzled as to why...

Thanks again, for the feed back and for Stan!

Well, if it matters at all I work for Mark, so I'm not coming from an orthogonal direction!

Firstly, the big issue here is sampling jointly or sampling conditionally.  Sampling conditionally is always going to bring in problems when the target measure has large correlations, i.e. any realistic problem.  HMC appears to use a Gibbs step to sample the auxiliary momenta, but technically that sampling is doing something else.  Ignoring the technical details, HMC is making a joint update of all the model parameters which is what is important.  

This leaves two issues: hierarchical modes are hard and discrete models are harder.

We review hierarchical models thoroughly in http://arxiv.org/abs/1312.0906, in particular we address why Euclidean HMC can run into problems.  A common approach that people take, starting with Radford Neal way back when he was first introducing HMC into stats, is to sample from the hyperparameters conditional on the local parameters,

but this inevitably leads to large autocorrelations because of the strong dependencies between the hyperparameters and the local parameters.  What you'd really want to do is sample from the marginal distribution of the hyperparameters, but computing that distribution is a tremendously difficult task.  Pseudomarginal methods are trying to solve this, but they are hard to make general and can run into their own pathologies.  As noted in the paper, Riemannian HMC solves many of these problems by construction, and it's coming soon to Stan.  At the moment your best bet is to find a good non-centered parameterization.

Everyone wants discrete parameters but everyone also needs to understand the underlying problems.  Without the nice structure of a continuous measure, it's very, very hard to construct an efficient sampler.  Moreover, if you split your parameters into continuous and discrete then you're doing conditional updates, which will be slow as discussed above.  We've seen lots of success when the discrete parameters can be marginalized out to give a fully continuous model amenable to HMC.  This includes discrete latent models, provided there aren't too many possible states.

You can certainly use Stan as a sampling API and stick it into a bigger sampler, but without addressing the fundamental problems with discrete variables and conditional updates you're only going to be able to go so far.

[david...@gmail.com]

Thanks for the responses Michael,

I certainly defer to you, Girolami and the rest of the Stan team on these issues.  From a quick skim of the paper it seems the argument is that centring is only possible under limited circumstances...  I agree it is a good paper, and will go back to it when I get the chance...  so far I have read Murray's and Filippone's paper - rather than yours.  Filippone's paper recommends interweaving using HMC for latent variables and plain old random walk Metropolis for hyperparameters.

Obviously conditional Gibbs sampling or single component Metropolis sampling is a compromise - but its a compromise that has served some of us quite well at times.

I would like to implement the model described here.  http://people.ee.duke.edu/~lcarin/PoissonProcess.pdf
It doesn't have to be using their sampler which involves both Metropolis Hastings and HMC steps - and involves trans-dimensional steps - but my first thought was to just use their sampler, and I thought Stan could help with the HMC part - but right now this appears to be not easy. 

OTOH If you think it can be completely done in Stan then I would be fascinated!

BTW for reference it seems that you have been over this issue just a few weeks back (sorry) see Krzysztof Sakrejda's post and the responses which gives some detail about why the above doesn't work...
https://groups.google.com/forum/#!topic/stan-users/HBewhdW1u88

[Michael Betancourt]

I certainly defer to you, Girolami and the rest of the Stan team on these issues.  From a quick skim of the paper it seems the argument is that centring is only possible under limited circumstances...  I agree it is a good paper, and will go back to it when I get the chance...  so far I have read Murray's and Filippone's paper - rather than yours.  Filippone's paper recommends interweaving using HMC for latent variables and plain old random walk Metropolis for hyperparameters.

Random walk on the hyperparameters isn't going to do anyone any good.  Certainly in the worst case HMC reduces to a random walk, so why not just use it in the first place?

Non-centered parameterizations work really well in practice, even when there are multiple levels and formally a few different ways of centering.  One of the nice things about

Stan is that it makes experimenting with different parameterizations straightforward.

 

Obviously conditional Gibbs sampling or single component Metropolis sampling is a compromise - but its a compromise that has served some of us quite well at times.

Of course any algorithm will work well in some circumstances, but we're (okay, maybe I should just speak for myself) interested in algorithms that work well with great generality,

and Gibbs is not going to be that algorithm, especially for modern applied problems.

 

I would like to implement the model described here.  http://people.ee.duke.edu/~lcarin/PoissonProcess.pdf
It doesn't have to be using their sampler which involves both Metropolis Hastings and HMC steps - and involves trans-dimensional steps - but my first thought was to just use their sampler, and I thought Stan could help with the HMC part - but right now this appears to be not easy. 

OTOH If you think it can be completely done in Stan then I would be fascinated!

BTW for reference it seems that you have been over this issue just a few weeks back (sorry) see Krzysztof Sakrejda's post and the responses which gives some detail about why the above doesn't work...
https://groups.google.com/forum/#!topic/stan-users/HBewhdW1u88

Trans-dimensional models are a different problem altogether, and I don't know of any algorithm that works particularly well outside of knowing a reasonable truncation point a priori and building a model of constant dimension.  Then again, we haven't spent a tremendous amount of time trying to build and fit these models.

Interestingly enough there are some possible generalizations to algorithms like HMC to these models, but the math is really complicated and somewhat obscure, so they won't be coming any time soon.

[david...@gmail.com]

Random walk on the hyperparameters isn't going to do anyone any good.  Certainly in the worst case HMC reduces to a random walk, so why not just use it in the first place?

I believe its the overhead of HMC.  On a per iteration basis I agree random walk Metropolis has no advantage - but random walk Metropolis will iterate faster...  but see the conclusion of Filippone rather than have me paraphrase...

 

 Of course any algorithm will work well in some circumstances, but we're (okay, maybe I should just speak for myself) interested in algorithms that work well with great generality,

and Gibbs is not going to be that algorithm, especially for modern applied problems.

 

I take your point completely - but would still like to be able to use Gibbs/single component Metropolis for problems that I don't know how to do anything else with - and for which there is some sort of track record of at least moderate performance.... 

I really look forward to modifications to rstan that make this possible.  Thanks again for your advice and pointers!

Cheers,
David

[Nigel Goodwin]

I'm not a stan user, but I have implemented NUTS together with other samplers, particularly DREAM (and RWM). The reason is two-fold - DREAM can handle, at least to some degree, discrete parameters and multi-modal situations. Secondly, when DREAM works well, it works well, and can be faster than NUTS. So I mix randomly DREAM and NUTS, in concept NUTS is used to generate a good population, and then DREAM is used to sample and difference from that population. This is all done randomly at each step and for each chain.

 

This also avoids one of the main pitfalls of DREAM, IMHO, that heuristic tricks have to be used to make it go outside of the hyperplane defined by the chains at the current step. The inclusion of NUTS means that those tricks are no longer needed.

 

In addition, I examine the acceptance rate for DREAM, and vary the proportion of time spent in DREAM according to how close the acceptance rate is to the theoretical value for normal distributions. I use dual averaging for the DREAM step (i.e. the scaled difference).

 

Remarkably, the dual averaging generated a step size almost identical to the theoretical value for a normal test problem. Well, I'm always surprised when theory meets practice.

 

So, long story short, I found value in mixing samplers. It's the same when you mix optimisation algorithms, they all have different good and bad behaviours, and a mixture is more robust.

[Bob Carpenter]

The technical problem right now in plugging samples back into
Stan is that there's not an easy way to pull the adapted mass matrix out
and then plug it back in. We'll be making that easier going forward,
probably as part of Stan 3, which is not on the near term horizon yet,
but we're working on design.

- Bob

[Bob Carpenter]

How do you deal with the interleaving when DREAM assumes an ensemble
and HMC/NUTS don't? And just to be clear here, DREAM isn't a random
walk algorithm in the same way as Metropolis and Gibbs are.

We're working on some ensemble samplers now, starting with the Goodman-Weare
versions, but we'll probably also implement differential evolution (which in my
understanding is like DREAM w/o the detailed-balance violating heuristics).

- Bob

[Nigel Goodwin]

Bob, I'm afraid I don't know exactly what you mean by 'ensemble' and 'interweaving'. But maybe this will help:

 

- I use multiple chains (8, for no particular reason), and parallelise at the likelihood function level, not the chain level. So I loop over steps, and for each step I loop over chains. I choose at random (with sampler weights) for each step/chain what sampler I am using, so for a given step there are multiple samplers, and of course the NUTS sampler takes a lot longer than the DREAM sampler, but that is OK because of the way I parallelise. It would be hopeless if I parallelised at the chain level.

 

- there are several different versions of DREAM. Some sample from other chains at the same step, others sample from a store of previous steps (DREAMZ). I use both methods. In pure DREAM, the problem with sampling from other chains is the one I mentioned, getting stuck in the same hyperplane, and they have various methods to overcome this, none of which I like. They are Gibbs-like, which is hopeless if the likelihood is not aligned with the axes. But by combining with NUTS, this problem is avoided. All the complexity of DREAM (choosing subsets of parameters, adding in a bit of random walk) can be avoided.

 

- I have a feeling, and it is no more than that, that DREAM sampling from other chains has a more satisfactory detailed balance proof, hence why I use multiple chains. But I may be wrong. the 'proofs' in DREAM are rather sketchy.

 

- whether the DREAM family, or my combined approach, strictly obey the necessary rules, is not something I am qualified to discuss. But in my tests they successfully reproduce the analytical results for complex high dimensional problems, which is more than good enough for my audience.

 

- speed and performance is critical in my application, it is an interactive tool for engineers which repeats MC many times, so I have to make some compromises to get the performance. Whether the compromises are real or imaginary I don't know (see above!), but without performance nobody would use it.

 

My video http://www.youtube.com/channel/UCQbRkgaevt0TSDzipvmtMKQ

 

shows, for EssRisk64DNormal, this combined approach for a 64 dimensional multivariate normal. In this case I have one chain, and DREAM samples from the store. Blue jumps are DREAM, green jumps are NUTS.

 

For the Haario problem (64 dimensional, each pair of parameters is a 2D Haario function, so it is much harder than the published function) DREAM works much less well, for obvious reasons. But my algorithm can recognise this, and adjust the amount of DREAM usage accordingly.

 

In concept, I guess you could say NUTS does the exploration, and other samplers fill in the gaps.

 

Also, when you combine samplers, you may need to look at modifying the algorithms for the individual samples, as I did for DREAM.  This is just the same when, for example, you combine evolutionary algorithm optimisers with local gradient optimisers - just plugging the two together is not a satisfactory solution, they have influence on each other, and opportunities also.

 

I really must look at Goodman-Weare in more detail one day.

 

Nigel

[Bob Carpenter]

What I meant is how do you combine HMC and ensemble methods.

An ensemble method uses K ensemble members, each being
a complete set of parameter values, whereas HMC works for a
single value of the parameters. The entire ensemble is what
makes up the Markov chain in the ensemble method --- the
individual ensemble members do not form a Markov chain by
themselves.

You can find the detailed balance proofs for the simplified
version in ter Braak's paper on differential evolution.

- Bob

[Krzysztof Sakrejda]

On Tuesday, July 22, 2014 12:41:40 PM UTC, david...@gmail.com wrote:

Thanks for the responses Michael,

I certainly defer to you, Girolami and the rest of the Stan team on these issues.  From a quick skim of the paper it seems the argument is that centring is only possible under limited circumstances...  I agree it is a good paper, and will go back to it when I get the chance...  so far I have read Murray's and Filippone's paper - rather than yours.  Filippone's paper recommends interweaving using HMC for latent variables and plain old random walk Metropolis for hyperparameters.

Obviously conditional Gibbs sampling or single component Metropolis sampling is a compromise - but its a compromise that has served some of us quite well at times.

I would like to implement the model described here.  http://people.ee.duke.edu/~lcarin/PoissonProcess.pdf
It doesn't have to be using their sampler which involves both Metropolis Hastings and HMC steps - and involves trans-dimensional steps - but my first thought was to just use their sampler, and I thought Stan could help with the HMC part - but right now this appears to be not easy. 

Thanks for sharing that, I was not aware of the sampling of event times by thinning.  If you are willing to take a step back from the full GP version and don't need to simulate events from the IHPP it looks like Stan can work just fine.  I've done a little proof of concept that this is reasonable in Stan by representing the intensity function as a spline with smoothness provided by a GP on the knot weights.  That gets you an arbitrarily fine intensity function and reduces the dimension of the GP to the number of knots (rather than the number of data points).  Once you estimate the spline you could even simulate from the poisson process by thinning as suggested by Adams et al.  I'm attaching some code so it's a little more concrete, not that I'm claiming this to be a _good_ model for the simulated data but Stan does pretty well, all things considered.
 

OTOH If you think it can be completely done in Stan then I would be fascinated!

BTW for reference it seems that you have been over this issue just a few weeks back (sorry) see Krzysztof Sakrejda's post and the responses which gives some detail about why the above doesn't work...
https://groups.google.com/forum/#!topic/stan-users/HBewhdW1u88

I appreciate it getting re-hashed in more detail.  In my case I think one of the solutions is not sticking to discretized time survival models quite so much.

Krzysztof
 

[Nigel Goodwin]

OK, I think I understand, the ensemble is the set of chains.

 

As I say, I trample all over the theory, in an attempt to get something which is useable by engineers, and much better than what they had before (RWM).

 

I think I can see your worry - but it will need some further thought to see if it clearly breaks the detailed balance proof. I have of course read most of the DREAM papers and tried to understand the detailed balance proofs. I still have an ensemble, but different members use different samplers at different steps.

 

DREAMZ does not, as far as I recall, rely on the ensemble for detailed balance proof, it uses a store of previous samples.

 

NUTS by itself is too slow for me and doesn't handle discrete cases, and my implementation has been in commercial use for over a year now, and cash flow is king. No doubt something better will emerge over time. For at least the next year my focus will be on dissemination rather than major changes to algorithms.

[david...@gmail.com]

Thanks for your code Krzysztof, I didn't expect this!

It doesn't run out of the box for me, I will try tweaking it to try to work out what you have done.

I hadn't thought about using splines, but it is not clear to me if doing so solves everything.  Maybe you can analytically integrate under a non-negative function of them which would change things....

For this model, the likelihood contains the following terms

1 ~ poisson(sigmoid(lambda))
0 ~ poisson(sigmoid(-latent))

lambda is the size of the data set, where latent is of variable size (which is the big problem for Stan I think).

The reason for the use of the 1/0 is that for a point process there is a small number of locations where a single event happens and everywhere else nothing happens - you never evaluate the poisson with a count greater than 1.

it looks to me like you are instead modelling

count[ii] ~ poisson(lambda(ii))

which appears at face value to be a different model.  Apologies, in advance if I have misunderstood...  Maybe we should talk more about this out of the thread...  my email is david...@gmail.com.

Cheers,
David

[Krzysztof Sakrejda]

On Thursday, July 24, 2014 7:51:35 AM UTC, david...@gmail.com wrote:

Thanks for your code Krzysztof, I didn't expect this!

You're welcome

It doesn't run out of the box for me, I will try tweaking it to try to work out what you have done.

I think there is one hard-coded path in there for loading the file and maybe something else that needs tweeking but if after that you're actually getting errors I'll recheck my files and upload more carefully.

I hadn't thought about using splines, but it is not clear to me if doing so solves everything.  Maybe you can analytically integrate under a non-negative function of them which would change things....

The nice thing about using a gaussian basis for the splines with weights bounded positive is that you can integrate analytically and that's what I did for this code.  That's why the expression for lambda is in terms of a sum of CDF's. I'm pretty sure this is straightforward using any pdf with a nice CDF for a basis.
 

For this model, the likelihood contains the following terms

1 ~ poisson(sigmoid(lambda))
0 ~ poisson(sigmoid(-latent))

lambda is the size of the data set, where latent is of variable size (which is the big problem for Stan I think).

ok.

The reason for the use of the 1/0 is that for a point process there is a small number of locations where a single event happens and everywhere else nothing happens - you never evaluate the poisson with a count greater than 1.

ok, that's what you'd get for a continuous time survival model too, I don't see why count can't be zero/1... 

it looks to me like you are instead modelling

count[ii] ~ poisson(lambda(ii))

which appears at face value to be a different model.  Apologies, in advance if I have misunderstood...  Maybe we should talk more about this out of the thread...  my email is david...@gmail.com.

I don't see why one can't be turned into the other.  I think when you code a continuous time survival model based on a poisson process you end up with something like you are doing but it should
still be fine with modification of the code.  That's on my to-do list so I will probably post both of these to the list to get comments in a separate thread (I guess you're right we're off-topic here)

Krzysztof
 

[Nigel Goodwin]

Looking again at the Braak 2006 Stat Comp paper,he says quote ".... and the distribution of e is symmetric........For the third condition it is required that any state can be reached with positive probability and this is guaranteed by the unbounded support of the distribution of e (which is N(0, b))" Now if e is a NUTS step (or any other valid sampler) rather than a RWM step, does this invalidate the proof? It isn't obvious to me, but I'm not an expert. Is the problem that the NUTS is not a symmetric distribution?

 

I somehow thought that if you have an ensemble, so long as each individual step conformed to the rules, the ensemble would conform to the rules, so it didn't matter if individual steps and members of the ensemble used different samplers.

 

I guess clarification of this would help those of us doing what this thread is about, using NUTS as a sampler within a bigger algorithm.

[Nigel Goodwin]

ps. Michael's paper on thermodynamic HMC looks very interesting. I have played around with tempering, and anything which helps with adaptive and automatic tuning is a good thing - that is why NUTS is such a breakthrough compared with previous HMC methods.

 

 

[Bob Carpenter]

Yes --- you can interleave steps from Markov chains
with the right stationarity properties (i.e., the
stationary distribution of the chain is the posterior
p(theta)). So you can alternate differential evolution
and RWM or HMC or NUTS and the combination still has
the right stationarity properties.

The detailed balance proof for NUTS is rather complicated,
but similar in many ways to the arguments involved for
slice sampling, which Radford Neal goes over in his paper
introducing slice sampling.

- Bob

On Jul 24, 2014, at 4:00 PM, Nigel Goodwin <goo...@essence-property.com> wrote:

> Looking again at the Braak 2006 Stat Comp paper,he says quote ".... and the distribution of e is symmetric........For the third condition it is required that any state can be reached with positive probability and this is guaranteed by the unbounded support of the distribution of e (which is N(0, b))" Now if e is a NUTS step (or any other valid sampler) rather than a RWM step, does this invalidate the proof? It isn't obvious to me, but I'm not an expert.

Everything should still be OK. The theory is that you need
to ensure that the stationary distribution of the Markov chain
is the posterior p(theta). If that holds, as it does for
differential evolution and NUTS or HMC (or RWM), then you can


- Bob

[david...@gmail.com]

If you have managed to find an analytic integral that removes the doubly intractable problem for even a grossly simplified model that would be quite remarkable...

Its not a problem to integrate under a standard GP (for some covariance functions e.g. gaussian), but once you apply a non-negative transform (log, square or sigmoid) it becomes intractable AFAICT and right now I don't see how using a spline will help... I can't spend too much time on this right now, but we should keep talking about this at a later date...

Cheers,
David

[Nigel Goodwin]

Thanks for the feedback, it seems I'm not obviously doing anything ridiculous, and alternating samplers is, in principle, OK.

 

Just like the no free lunch theorem (which I've seen used to justify all kinds of nonsense) in optimisation, different samplers may be suitable for different problems, and then the trick is to automatically adjust the proportion of the time spent in different samplers according to the problem. This was the main trick I used with DREAM, by examining the step size from dual averaging v. the theoretical step size for normal distributions. In practice, it means that if the hessian is fairly stable, DREAM dominates, if not NUTS dominates. It is easy to see why DREAM is hopeless with severe bananas..

[Bob Carpenter]

You just have to be careful to do all this balancing
during warmup, throw away the warmup samples, and then
used fixed parameters during actual sampling.

If you adjust ratios after warmup, then you break
detailed balance and almost certainly the
break the stationary distribution for MCMC, except
under circumstances where the adjustment tails off
fast enough in magnitude or the adjustments themselves
are reversible.

- Bob

[Nigel Goodwin]

That's what I do.