Toonces The Driving Model: Quasi-discrete model does converge, just not very well!

[Phil Price]

This is in some ways a follow-up to a post I had a month ago, but it's different enough that I think it justifies its own post. 

I have created a super simple example of what I am trying to do. I have a time series that has a single step in it. Basically I want to find the start and end of the step, and the magnitude. In a real-world case there would be noise added, plus a bunch of other things that are far more complicated (see earlier post for the complications if you interested) but for this test problem it is really easy. I have attached a file for one-stop shopping: source the file in R (and wait a while) and it should just work.

As written, the data vector y has a series of 1 for a while; then they become 1.6 for a while; then down to 1 again. 

This is seemingly a discrete problem: there is some integer index istart which is the first place y[istart] becomes 1.6, and another index iend which is the last place where y[iend] = 1.6 before dropping down to 1. 

What happens if I define continuous variables start and end that define the range? In my earlier problem this did not work. The issue is that if there's no difference in the model predictions when start = 10.1 vs start = 10.9 (for example), the likelihood function doesn't vary continuously and the sampler gets stuck. In either case, the prediction takes one value up to i=10 and another prediction from i=11 onward. 

I fixed this by adding something to the prediction _before_ the start time. If the start time is 10.9 I add very little to the point at i=10; if the start time is 10.1 I add a lot to the point at i=10.  Now, if the real start time is i=11 and the sampler is at 10.1, the likelihood will improve monotonically as the sampler moves towards 11.  Genius! (I do the same thing at the other end. By the way, I implement this as a start time and a duration, not a start time and an end time). 

The good news is: it works.

The bad news is: it works very very slowly. 100 or even 1000 iterations won't do it -- indeed, I initially thought it wasn't working at all, until I looked at the chains and saw very slow progress towards the right answer. At 100000 iterations (of which the first half are warm-up) it is clearly making good progress, but still hasn't gotten down to a narrow interval around the correct answer. 

After posting this, I will let it do 200000 and see where that gets me. I will post a follow-up here.

If convergence takes hundreds of thousands of iterations even with this simple problem, there's no way I can use this approach for my real problem, which will have dozens of steps at unknown times, as well as other complications. I'm hoping there's an easy improvement.

I notice that the sampler takes really small steps. It will try a start of 20.0133, and then the next iteration it will be at 20.0131 or something. Obviously it will take a long time to reach the right answer of 5, even if it moves the right direction the whole time. 

I'm aware that there are some stan parameters like adapt_delta and adapt_kappa and adapt_window and so on, that might matter here. But I haven't been able to find the documentation that explains these, and even if setting these parameters is the solution to my problem I could use some wisdom in how to figure out which ones to monkey with and in which directions, and, equally important, how to figure it out for a more complicated problem like the one I am working towards.

As always, advice is greatly appreciated. 

Phil Price

[Phil Price]

I found a bug in my model. Fixing it doesn't help, but its existence does make me question whether there are others. Please don't spend time looking at the code until I post again with code I'm sure is right.

Phil

[Bob Carpenter]

Have you looked at the change-point model in the manual? It's
the first example in the latent discrete parameters chapter.
That'll be the most efficient way to code this kind of thing.
But it'll still be cubic in the number of points if you have
two change points, even if the second one is by construction returning
to the original value.

Andrew always complains about these examples and says they're
unrealistic --- that we should just be using an ordinary time series
kind of model because any point can change a lot or a little. Maybe
a mixture noise model on the time series to allow local and big jumps.
Then let the data sort it out. That's a different model, but it
won't be combinatorially terrible.

Coding this up with a continuous parameter to approximate a discrete
parameter isn't likely to be efficient, even if you can get it to work.
I see you tried to hack it so derivatives went the right way, but
I have no idea how to set that up in general. I think Andrew and I are sending you conflicting advice here.

- Bob

> On Dec 29, 2016, at 2:34 PM, Phil Price <pnp...@creekcats.com> wrote:
>
>
> I found a bug in my model. Fixing it doesn't help, but its existence does make me question whether there are others. Please don't spend time looking at the code until I post again with code I'm sure is right.
>
> Phil
>
>

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Phil Price]

On Thursday, December 29, 2016 at 2:32:29 PM UTC-8, Bob Carpenter wrote:

Have you looked at the change-point model in the manual?  It's
the first example in the latent discrete parameters chapter.
That'll be the most efficient way to code this kind of thing.
But it'll still be cubic in the number of points if you have
two change points, even if the second one is by construction returning
to the original value.  

I haven't looked at it, but will.

 

Andrew always complains about these examples and says they're
unrealistic --- that we should just be using an ordinary time series
kind of model because any point can change a lot or a little.  Maybe
a mixture noise model on the time series to allow local and big jumps.
Then let the data sort it out.  That's a different model, but it
won't be combinatorially terrible.

This is my approach in my original model -- the one I showed the Stan group -- and I am 90% happy with that model. But it doesn't really handle big step changes well (as you know but other readers don't). Uh...this won't be useful to anyone unless I give a lot of background, which I won't do here. I'll just say that since one jump of 2 standard deviations is much less likely than two consecutive jumps of 1 standard deviation, my previous implementation would tend to smear step changes over several time intervals, and that is undesirable for my application. So I'm trying to implement a model that allows real step changes. Depending on how well I can get it to work, I may completely rewrite my model so that it _only_ uses step changes, allowing different durations and lengths...perhaps many many changes... or maybe I'll do a hybrid approach, where I modify my orignal model to allow N step changes in addition to the way it usually handles stuff.

Coding this up with a continuous parameter to approximate a discrete
parameter isn't likely to be efficient, even if you can get it to work.
I see you tried to hack it so derivatives went the right way, but
I have no idea how to set that up in general. I think Andrew and I are sending you conflicting advice here.

I ended up with a great illustration of Andrew's "folk theorem" that if a model is very slow to converge, you probably have a coding error. I had an off-by-one error in the way I was handling my trick to make the derivatives go the right way. Once I fixed that, it works great, super fast convergence on my toy problems. I will build more and more complicated toy problems until I am ready to try my real one, but for now it is promising. 

Thanks, Stan community, for being a rubber duck.  Thanks, Bob, for your suggestions. 

I am attaching a fixed version of my toy problem, in case anyone wants to take a look.

By the way, Bob (or others on the Stan team) I'd like to put in a plug for adding the functions floor(), ceiling(), and maybe round() to the language. I am mimicking these with a while loop that counts through integers until it exceeds the number I'm interested in -- and it was an embarrassing mistake in doing this that led to my off-by-one problem, but of course that's no one's fault but mine. But it's still a pain to have to program this everywhere I need it, and it hurts readability, and...well, basically there's a good reason pretty much every other language lets you find the integer part of a number.  

Also, although it's not a big deal to count to 90 when trying to find the greatest integer less than 90.3, it is certainly inefficient...and what if I am working with, say, a year of 1-second data (not at all an unrealistic possibility...indeed, I may soon be working with years worth of 1/120-second data!). I'm not going to count to 27 million to find the integer part of 27,000,000.63.  There must be a super-fast way to do this... find a range and then use bisection, perhaps? Maybe there's something really simple that should be obvious to me but isn't. But why make me jump through any of these hoops at all? Please, put real --> int conversion into the language.

Phil

 

[Andrew Gelman]

> On Dec 30, 2016, at 2:02 PM, Phil Price <pnp...@creekcats.com> wrote:
>
> Once I fixed that, it works great, super fast convergence on my toy problems. I will build more and more complicated toy problems until I am ready to try my real one, but for now it is promising.

Great!

> By the way, Bob (or others on the Stan team) I'd like to put in a plug for adding the functions floor(), ceiling(), and maybe round() to the language.

Huh? Stan has these functions! See page 420-421 of the Stan 2.14.0 manual! (And the functions are not new to version 2.14.)

[Bob Carpenter]

> On Dec 30, 2016, at 2:02 PM, Phil Price <pnp...@creekcats.com> wrote:
>
>
> By the way, Bob (or others on the Stan team) I'd like to put in a plug for adding the functions floor(), ceiling(), and maybe round() to the language. I am mimicking these with a while loop that counts through integers until it exceeds the number I'm interested in -- and it was an embarrassing mistake in doing this that led to my off-by-one problem, but of course that's no one's fault but mine. But it's still a pain to have to program this everywhere I need it, and it hurts readability, and...well, basically there's a good reason pretty much every other language lets you find the integer part of a number.

This is an ongoing request we receive all the time. There
are warnings in the manual around places where we do allow it.

We do allow you to get the floating point integer part
of a number, but that's not useful for indexing. And even
that we shouldn't do, as the derivative's a step function
whereas HMC is assuming a continuously differentiable log density.

That's why we didn't add the versions as is. We may get rid
of the ability to do what you're trying to do in the future.

> Also, although it's not a big deal to count to 90 when trying to find the greatest integer less than 90.3, it is certainly inefficient...and what if I am working with, say, a year of 1-second data (not at all an unrealistic possibility...indeed, I may soon be working with years worth of 1/120-second data!). I'm not going to count to 27 million to find the integer part of 27,000,000.63. There must be a super-fast way to do this... find a range and then use bisection, perhaps? Maybe there's something really simple that should be obvious to me but isn't. But why make me jump through any of these hoops at all? Please, put real --> int conversion into the language.

Binary search will do it in log(n) time. We don't want to
make it hard and provide hoops, we want to discourage it altogether
in favor of continuous models (or marginalized discrete models).

Of course, everyone who wants to code this kind of thing just tells
us we should make it easier, as you did.

- Bob

[Bob Carpenter]

> On Dec 30, 2016, at 5:19 PM, Andrew Gelman <Gel...@stat.columbia.edu> wrote:
> ...

>> By the way, Bob (or others on the Stan team) I'd like to put in a plug for adding the functions floor(), ceiling(), and maybe round() to the language.
>
> Huh? Stan has these functions! See page 420-421 of the Stan 2.14.0 manual! (And the functions are not new to version 2.14.)

They return real values, not integers, so you can't use
them as indexes. That's intentional. Each one is linked
to a warning because it doesn't have continuous derivatives.
Which means it will break HMC unless used very carefully (as
in some kinds of interpolations).

- Bob