Extremely slow (or no) convergence of model with discrete steps...how to improve?

[Phil Price]

I am working with time series data. The actual application (to electricity consumption in large buildings) involves time series that have a regular daily pattern (low at night, high during the day) that is sometimes changed; for instance, lights might be scheduled to turn on at 6 a.m. daily until,  they are changed to turn on at 5 a.m. instead for a few weeks, then set back to 6 a.m.

The model I really want to fit is quite a bit more complicated than this, but let's start here.

I have tried an implementation of this that "should work" but doesn't. The problem seems to be with convergence rather than with the model per se. 

I have attached a file called "CreateData.R". Sourcing this file in R will create a time series y, as well as some ancillary parameter values (such as the number of data points per day, and their time-of-day indices).  This fake data set is very unrealistic: it is flat at y=1 every hour of every day for three weeks, then goes up to y=2 for the next three weeks, then back to y=1. Then some iid normal noise is added.

The model writes 

ypred[ t ] = betaTOD[ iWeek[t], iTOD[t]] 

where iWeek[t] is the week number ( = 1 for every hour of the first week, 2 for the second week, etc.)

and iTOD[t] is the time of day for time t.

the betaTOD values are themselves written as a sum:

betaTOD[jWeek, jTOD] = betaTOD[1, jTOD] + phiTOD[jWeek, jTOD]

That is, betaTOD in a given week is whatever it was in week 1 at that time of day, plus a change.

The change is itself a sum over "changes", where each change has a start time, a duration, and a magnitude delta:

 

 phiTOD[week, jTOD] = phiTOD[week, jTOD] + deltaTOD[iChangeTOD, jTOD];

 

 In the test example given, there is only one change for each hour of the day as determined  with nChangesTOD = 1. For the test data, this is the right number.

SimpleModel.Stan.txt contains a Stan model that fits the data. You can run it with a command like:

res = stan("./SimpleModel.Stan.txt", iter=400, chains=3, data=c("cooked", "nN", "nChangesTOD", "nTOD", "nWeeks", "iTOD", "iWeek", "sigmaMax", "y"))

The parameter "cooked" chooses between two models:

If cooked = 1, very tight prior distributions are set, for instance:

 if ( cooked==1 ) {

        weekStartTOD[iChangeTOD, jTOD] ~ normal(3.1,0.01);

        tauTOD[iChangeTOD, jTOD] ~ normal(3.1, 0.01);

        deltaTOD[iChangeTOD, jTOD] ~ normal(1.0, 0.1);

}

If cooked is not 1, then more real-world priors are used.

So here's the problem: if you use cooked = 0, you don't get the right answer...or, rather, maybe you can but you would need to set the number of iterations to some very very large number. 

But if you set cooked=1, you do get the right answer. (For this test problem, with only a small amount of noise added, the 'right answer' is quite tightly defiined, with distant parameter values being very unlikely.

What's happening is this: the "start time" and the "duration" are effectively quantized with units of weeks. They are defined as continuous variables, but whether the value of the start time is 2.1 or 2.9 or 2.99 makes no difference: all of these define the start time to be in "week 2". This makes it easy for the system to get trapped in a local likelihood maximum and unable to find the global one. 

Suggestions would be greatly appreciated!

Phil Price

[Bob Carpenter]

On Nov 26, 2016, at 9:03 PM, Phil Price <pnp...@creekcats.com> wrote:
>
>
> I am working with time series data. The actual application (to electricity consumption in large buildings) involves time series that have a regular daily pattern (low at night, high during the day) that is sometimes changed; for instance, lights might be scheduled to turn on at 6 a.m. daily until, they are changed to turn on at 5 a.m. instead for a few weeks, then set back to 6 a.m.
>
> The model I really want to fit is quite a bit more complicated than this, but let's start here.
>
> I have tried an implementation of this that "should work" but doesn't. The problem seems to be with convergence rather than with the model per se.

...

> What's happening is this: the "start time" and the "duration" are effectively quantized with units of weeks. They are defined as continuous variables, but whether the value of the start time is 2.1 or 2.9 or 2.99 makes no difference: all of these define the start time to be in "week 2". This makes it easy for the system to get trapped in a local likelihood maximum and unable to find the global one.

In general, that approach shouldn't work. The assumption built into
Stan's algorithms is that the log density function is continuously
differentiable. By hacking in discrete variables the way you do, you break
that property.

We see a lot of people try to do things like sample a Bernoulli(theta) by

parameters {
real<lower=0, upper=1> u1;

...

transformed parameters {
int z1 = u1 < theta;


This won't work. I think the impression people get is that Stan is
somehow sampling u1 first uniform(0, 1) and then assigning z1. But
that's not how it works. Instead, the algorithm feeds in values of
the parameters to evaluate, and Stan follows the gradients. Here there
are no gradients to follow.

The only thing that'll work well in Stan is marginalizing out the
discrete parameter. As we discussed when you're in town, this can be
combinatorially prohibitive. If it is, then you're just out of luck with
the model as formulated in Stan. Andrew or others may have ways to reformulate
to get the same effect, but I don't have any ideas.

- Bob

[Phil Price]

> What's happening is this: the "start time" and the "duration" are effectively quantized with units of weeks. They are defined as continuous variables, but whether the value of the start time is 2.1 or 2.9 or 2.99 makes no difference: all of these define the start time to be in "week 2". This makes it easy for the system to get trapped in a local likelihood maximum and unable to find the global one.

In general, that approach shouldn't work.  The assumption built into
Stan's algorithms is that the log density function is continuously
differentiable.  By hacking in discrete variables the way you do, you break
that property.  

I'm not surprised, indeed this is sort of what I expected...except that at the Stan meeting I attended a few weeks ago, I thought this was the model the Stan experts said I should fit.  

I thought it seemed possible that the prior distributions could provide a sort of tug in one direction or another: even though the predictions are exactly the same for any start times that have the same integer week number, the likelihood will not be the same. But you're right that it will be discontinuous across the boundary, so I guess that is the problem.

Did I misunderstand what Andrew and Michael were suggesting? 

 

The only thing that'll work well in Stan is marginalizing out the
discrete parameter.  As we discussed when you're in town, this can be
combinatorially prohibitive.  If it is, then you're just out of luck with
the model as formulated in Stan.  Andrew or others may have ways to reformulate
to get the same effect, but I don't have any ideas.

If I have to try every combination of start time and duration, that would be combinatorially prohibitive. With 100 weeks of data, there are 100 possible start times, and, uh...let's see, duration has to be less than or equal to the number of remaining weeks, so... that's (100*99) + (99*98) + (98*97) + ... different possibilities, which is 1 million. And each time of day could have a different start time and duration, so if I do hours that's 24 million (and in fact it would be more convenient to work with 15-minute data, so that's 100 million)...and this all assumes just one "impulse" per time of day, over a two-year period. A more common real-world case would be 10 or more.  10 would not be 10x as much computation, since we could set a constraint that the next impulse has to be after the first one has ended, but I would still be in the hundreds of millions of combinations. This is not an appealing approach. 

However, it seems like there should be a way to greatly reduce the amount of computation because it shouldn't be necessary to test every possible start time and duration:  Although p(y | start time, duration, magnitude ) is not a _continuous_ function of start time and duration, it is monotonic. If I were writing a search algorithm rather than a statistical model, I could find the week that is most different from the average and set the start time at the start of that week and the duration to 1 week. Then try increasing the duration to 2 weeks. Keep increasing the duration until things stop getting better. Then try stepping the start time back one week while increasing the duration one week (so the end time stays the same). Keep doing this until things stop getting better. I'm done. Then, to do another impulse, do the same thing with the residuals from the first one. And so on. Repeat until I'm happy with how well things fit. 

Perhaps there are other approaches too. 

Hmm. I would love to be able to do all of this in Stan, but perhaps I will need to try some kind of ping pong approach: use one system to come up with the first set of impulse start times and durations, then have Stan fit all of the continuous parameters, then pop back to system 1 to come up with the next impulse times, etc. Ugh.  

Thanks for your help, Bob. 

 

[Bob Carpenter]

> On Nov 29, 2016, at 6:58 PM, Phil Price <pnp...@creekcats.com> wrote:
>
>
>
> > What's happening is this: the "start time" and the "duration" are effectively quantized with units of weeks. They are defined as continuous variables, but whether the value of the start time is 2.1 or 2.9 or 2.99 makes no difference: all of these define the start time to be in "week 2". This makes it easy for the system to get trapped in a local likelihood maximum and unable to find the global one.
>
> In general, that approach shouldn't work. The assumption built into
> Stan's algorithms is that the log density function is continuously
> differentiable. By hacking in discrete variables the way you do, you break
> that property.
>
>
> I'm not surprised, indeed this is sort of what I expected...except that at the Stan meeting I attended a few weeks ago, I thought this was the model the Stan experts said I should fit.

> I thought it seemed possible that the prior distributions could provide a sort of tug in one direction or another: even though the predictions are exactly the same for any start times that have the same integer week number, the likelihood will not be the same. But you're right that it will be discontinuous across the boundary, so I guess that is the problem.
>
> Did I misunderstand what Andrew and Michael were suggesting?

Perhaps --- I'd ask them. I have a hard time imagining Michael
would've recommended this kind of discontinuous approach since
he's been arguing we take out the means to do it from the language
altogether (conditioning on parameters).

> The only thing that'll work well in Stan is marginalizing out the
> discrete parameter. As we discussed when you're in town, this can be
> combinatorially prohibitive. If it is, then you're just out of luck with
> the model as formulated in Stan. Andrew or others may have ways to reformulate
> to get the same effect, but I don't have any ideas.
>
>
>
> If I have to try every combination of start time and duration, that would be combinatorially prohibitive. With 100 weeks of data, there are 100 possible start times, and, uh...let's see, duration has to be less than or equal to the number of remaining weeks, so... that's (100*99) + (99*98) + (98*97) + ... different possibilities, which is 1 million.

That's what I was trying to say over lunch. Yes, that is
a problem.

> And each time of day could have a different start time and duration, so if I do hours that's 24 million (and in fact it would be more convenient to work with 15-minute data, so that's 100 million)...and this all assumes just one "impulse" per time of day, over a two-year period. A more common real-world case would be 10 or more. 10 would not be 10x as much computation, since we could set a constraint that the next impulse has to be after the first one has ended, but I would still be in the hundreds of millions of combinations. This is not an appealing approach.

Agreed.

> However, it seems like there should be a way to greatly reduce the amount of computation because it shouldn't be necessary to test every possible start time and duration: Although p(y | start time, duration, magnitude ) is not a _continuous_ function of start time and duration, it is monotonic. If I were writing a search algorithm rather than a statistical model, I could find the week that is most different from the average and set the start time at the start of that week and the duration to 1 week. Then try increasing the duration to 2 weeks. Keep increasing the duration until things stop getting better. Then try stepping the start time back one week while increasing the duration one week (so the end time stays the same). Keep doing this until things stop getting better. I'm done. Then, to do another impulse, do the same thing with the residuals from the first one. And so on. Repeat until I'm happy with how well things fit.

You could try to do something like Metropolis with a custom proposal
for the discrete parameters. We haven't played with this kind of
thing and don't know how changing the discrete parameters will
affect the sampling for the continuous parameters (the log density of
the continous parameters and hence the potential energy for HMC,
changes with the discrete parameters).

> Perhaps there are other approaches too.
>
> Hmm. I would love to be able to do all of this in Stan, but perhaps I will need to try some kind of ping pong approach: use one system to come up with the first set of impulse start times and durations, then have Stan fit all of the continuous parameters, then pop back to system 1 to come up with the next impulse times, etc. Ugh.

It's called "blocked" sampling and that's what PyMC3 does.
And NIMBLE's also set up to do that but I don't think they
give you HMC for the continuous parameters.

- Bob

[Michael Betancourt]

Perhaps --- I'd ask them.  I have a hard time imagining Michael
would've recommended this kind of discontinuous approach since
he's been arguing we take out the means to do it from the language
altogether (conditioning on parameters).

I don't remember ever talking about discrete parameters or associated mixtures.  We were talking about time series with multiple time components identified by priors.  For example,

f(t) = A1 sin(t / T1) + A2 sin(t / T2) + A3 sin(t / T3);

T1 ~ normal(0, 20 minutes);

T2 ~ normal(1 hour, 5 minutes);

T3 ~ normal(1 day, 2 hours);

In particular you were talking about boxcar-like effects that could turn on or off for different periods, which could be modeled with a continuous start time (perhaps guided towards certain times of the day by priors) and a continuous duration (identified into time scales by priors like above).  

But at no point were we talking about trying to turn the components on or off -- you just assume that they always there and if the data don't need them then they'll be pushed towards zero and have minimum effect on the fit.

 

[Andrew Gelman]

Hmm, there are a couple of ways I can thinking of hacking this (or maybe "solving" rather than "hacking" if we can consider these changes as representing a true underlying model):

1.  Continuous-time parameter, continuous-time data:  Allow these weekly effects to start and end midweek.  For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.  With real data, if things actually only start at the beginning of each week, you might still be able to fit this model because the slop and uncertainty would allow a continuous fit.  Sure, it wouldn't be quite correct, but you already have enough noise that you're not always identifying the start and end _weeks_, so would it really be so horrible to allow continuous times?

2.  Continuous-time parameter, discrete-time data:  Again allow a continuous time for the start and end parameters.  But now match them to discrete weekly start and end via an S-curve.  For example, let's label the discrete points as integers.  Then if the continuous start parameter is x, the discrete starting point could be the integer closest to x + epsilon, where epsilon is a distribution such as a logistic or normal with center 0 and scale sigma, and sigma has some value such as 0.5 or 1:  this corresponds to an ordered logit or probit.  The idea is that sigma will be small enough that the start time parameter roughly localizes the actual start time, but large enough that the transition from one integer to the next is smooth, so there's no near-discontinuity that would jar the HMC.  You could implement this model with latent epsilon but I think it would probably work better to use the normal and logistic cdfs.  (In theory these go off to infinity but I'd just truncate to the nearest 3 or 4 discrete possibilities for the likelihood.)  Once you've fit the model, you can then, in the generated quantities block, use simple Bayes calculation to get inference for the epsilons and thus identify the actual starting and end points more precisely.

Of these, I think I prefer #1.  But I'd be happy with either one if it works!

See you

A 

--
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]

I don't remember ever talking about discrete parameters or associated mixtures.  We were talking about time series with multiple time components identified by priors.  For example,

f(t) = A1 sin(t / T1) + A2 sin(t / T2) + A3 sin(t / T3);

T1 ~ normal(0, 20 minutes);

T2 ~ normal(1 hour, 5 minutes);

T3 ~ normal(1 day, 2 hours);

In particular you were talking about boxcar-like effects that could turn on or off for different periods, which could be modeled with a continuous start time (perhaps guided towards certain times of the day by priors) and a continuous duration (identified into time scales by priors like above).  

This example with the sin functions is far enough divorced from my case that I am having trouble mapping it onto my problem. In my problem, there are patterns that recur every 24 hours exactly, starting on some start day and ending on some other day. They are not necessarily continuous functions with period 24 hours, although they can be. For instance, suppose there is exterior lighting that turns on every day at 6pm and turns off every day at 7am. I can already capture that recurring pattern. But suppose, starting November 1 and ending the next March 1, the lighting turns on every day at 5pm instead. Let's say the lighting is 1 kW. As a deviation from the original situation, this shows up as a new spike in the hourly data point at 5 pm. Everything else remains unchanged.  

There are different ways of conceptualizing or modeling that case, which looks like a square wave of off-on-off-on (or, if you like, a chain of boxcars) for which the length increases by an hour from November 1 to March 1. (Of course, in the real-world application, this square wave will be superimposed on the many other electric loads in the building). As described in the paragraph above, I picture it as a recurring square wave that doesn't vary in length, plus the one-hour pulse every day from Nov-Mar. But I'm happy to consider other ways of expressing it. However: sometimes there really is something that happens only as a single data point each day, starting on one day and ending on another: Every day from 1am-2am the industrial ovens are turned to full blast to clean them, things like that. So whatever the parameterization, it does have to be able to handle spiky behavior that starts one week and ends another.

My examples above assume a period of 24 hours. There can also be changes with a period of 1 week (e.g. something new happens every Monday at 8 am or whatever). My desired model includes both. 

I don't see how the paradigm you are suggesting would allow me to model this. Unlike your example, the periods are not unknowns to be estimated from the data. I know the periods, what I don't know is the start and stop times: when does the boxcar start and how long is it? My current (non-working) model makes the start time and duration continuous variables, but this doesn't fit the application: what does it mean to say that starting (week 4, Tuesday, 3pm) we add 1 kW daily at 5pm? It doesn't mean anything, and therein lies the problem.

See below for why I don't just let every start and stop time be in the mix, with only the magnitudes unknown.

But at no point were we talking about trying to turn the components on or off -- you just assume that they always there and if the data don't need them then they'll be pushed towards zero and have minimum effect on the fit.

The last paragraph (immediately above)  is how I started to implement the impulse model, immediately after our discussion, and I may return to it: I assumed that every day at every hour was the start of a new impulse of unknown duration and magnitude. So, for instance, a new impulse starts on week 3 at 4 pm that applies every day at 4pm until the end of the impulse. This will not work, however, because the duration can only be an integer number of weeks so I run into the same issue of non-differentiability that is plaguing my current implementation. There might also be an issue with selecting priors so as to capture what I want without overfitting, but perhaps I'd be able to resolve that.

Another option, at least conceptually, is to include every start time and duration so that only the magnitudes are unknown: Have an impulse that occurs every day at 3pm in week 1 with a duration of 1 week; another one that occurs every day at 3pm in week 1 with a duration of 2 weeks; and so on for every duration and every week and every time of day. There are some major identifiability problems if I try it this way, but perhaps they could be solved with a judicious choice of priors. The bigger problem is simply the enormous number of parameters, which would be well into the millions for a year-long dataset. 

So...sorry for being a dope, but I'm still not getting it.

[Phil Price]

On Thursday, December 1, 2016 at 10:33:16 AM UTC-8, Andrew Gelman wrote:

Hmm, there are a couple of ways I can thinking of hacking this (or maybe "solving" rather than "hacking" if we can consider these changes as representing a true underlying model):

1.  Continuous-time parameter, continuous-time data:  Allow these weekly effects to start and end midweek.  For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.  With real data, if things actually only start at the beginning of each week, you might still be able to fit this model because the slop and uncertainty would allow a continuous fit.  Sure, it wouldn't be quite correct, but you already have enough noise that you're not always identifying the start and end _weeks_, so would it really be so horrible to allow continuous times?

Yes, I've thought about this, at least in the context of the daily changes. Indeed, the "time of day" effect really should be able to change at any day, and not just change at the start of a week, so it would be good if it can change any day, not just any week. But this doesn't solve the problem. See below for a case that is super easy to see.

Consider the "time of week" effect -- something that happens, say, every Monday at 3 pm. It really does not make any sense to say that the "Mondays at 3pm" effect changed starting on Tuesday at 8:14.  When I say "it doesn't make any sense", what I mean is that my predictions will all be exactly the same whether the "Mondays at 3pm" effect changes on Tuesday at 8:14 or on Friday at 11:45.  So I end up with the discontinous derivatives that are a problem with my current implementation too.

This happens with the "time of day" effect too. Suppose the building gets set, so that starting July 1 the air conditioning comes on at 7am daily instead of 6am. That is going to have no discernible effect outside the range 6am-8am or so.  Whether we say that the change happened on July 1 at 5am or on June 30 at 2pm will have no effect on the prediction and thus will lead to a discontinous likelihood. 

2.  Continuous-time parameter, discrete-time data:  Again allow a continuous time for the start and end parameters.  But now match them to discrete weekly start and end via an S-curve.  For example, let's label the discrete points as integers.  Then if the continuous start parameter is x, the discrete starting point could be the integer closest to x + epsilon, where epsilon is a distribution such as a logistic or normal with center 0 and scale sigma, and sigma has some value such as 0.5 or 1:  this corresponds to an ordered logit or probit.  The idea is that sigma will be small enough that the start time parameter roughly localizes the actual start time, but large enough that the transition from one integer to the next is smooth, so there's no near-discontinuity that would jar the HMC.  You could implement this model with latent epsilon but I think it would probably work better to use the normal and logistic cdfs.  (In theory these go off to infinity but I'd just truncate to the nearest 3 or 4 discrete possibilities for the likelihood.)  Once you've fit the model, you can then, in the generated quantities block, use simple Bayes calculation to get inference for the epsilons and thus identify the actual starting and end points more precisely.

Hmm...you're thinking of sort of tricking the solver by working at the limits of its numerical precision? If its step size can never go below 1e-12, we can make it so 1e-12 corresponds to a change of a whole integer? This seems incredibly kliudgy but would be kind of fun if it works. Can it possibly work, though? I would say that if the answer is "yes" then you guys should make discrete models "legal" in Stan. 

My gut feeling is that this will not work. Bob, what do you think?

Phil

 

[Michael Betancourt]

You'll have to model a mixture of duration times (1 unit of time, 2 units of time, 3 units of time...) for each possible effect.  Because the effects are independent you can do the marginalizations one at a time for a complexity of around n_effects * n_time_units.  In order to ensure that you don't get a multi modality you'll have to ensure that the effects don't overlap, but that will only be so successful if you have sparse/weak data (was the measurement low due to expected variability or did the effect turn off and a new one start? etc).

But honestly if you're talking about modeling day effects over the course of an entire year you're going to strain the possible computation anyways.  The only way this is going to be feasible is if you add more structure to the model (i.e. not just a collection of effects but a few AC effects with specific structure that differentiates them from other effects, etc).  

--

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+unsubscribe@googlegroups.com.

[Andrew Gelman]

I don't see why a mixture would be needed.  Why not just model it all directly?

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Phil Price]

On Thursday, December 1, 2016 at 5:57:40 PM UTC-8, Andrew Gelman wrote:

I don't see why a mixture would be needed.  Why not just model it all directly?

Because Stan won't let you evaluate duration = 1, duration = 2, duration =3, ... directly.   The only way you can do it is to have a mixture of duration=1, duration=2, duration=3, ..., and estimate the mixture coefficients.  Ditto for the start days or start weeks. 

With 24 hours per day, 50 or 100 weeks, 5 or 7 days per week (I can probably ignore weekends and just work with weekdays), and all possible durations from 1 week to N weeks, and allowing at least a few impulses for any given time of day...it's just way too much if I try to do a giant mixture model. Which is really too bad because in most cases I'm interested in, some small number of modes are going to account for the lion's share of the likelihood and I would be happy to completely ignore the others.  If someting turns on every day at 7pm from November through March, the mode that has it turning on at 7pm from August through May is going to have essentially 0 weight. 

I think I'm just SOL. Perhaps my best bet is, as I mentioned before, some sort of ping pong model: Fit a Stan model with no impulses; examine the residuals (outside Stan) to come up with candidate impulse start times and durations for each time of day; now fit a model that includes those (with fixed values for the start times and durations); rinse; repeat.  

I encourage the Stan team to consider extending the language to allow me to do all of this within Stan. I mean, obviously that isn't going to happen in time for me to use it for this specific problem, I'm not expecting that. But if I have to play ping pong, I'd rather do it all within one language, and to have access to the Stan data structures, rather than have Stan finish and output a bunch of stuff on which I have to write R functions (or Python or whatever) to extract the data I need. Perhaps this is a totally stupid suggestion, you guys are doing a great job with the language and clearly know your business so if you've thought about this and concluded it is untenable then I assume you're right.  But if it _is_ possible to extend the language to allow this, count me as one vote saying it's a good idea.

[Bob Carpenter]

I'm not sure what it is you want Stan (RStan?) to implement.
What you're trying to do sounds very particular to your specific
model. We don't want to try to take over R's job on the data
munging side. For example, I'm not sure what we'd provide
to "examine the residuals" that wouldn't be very model specific.
RStanArm can provide a lot of this kind of functionality because
it constrains the model classes and thus knows what everything
is.

- Bob

[Andrew Gelman]

Why not just let duration be a continuous parameter?

[Phillip Price]

On Dec 2, 2016, at 1:26 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Why not just let duration be a continuous parameter?

Duration can be a continuous parameter, and in my implementation (I posted the code), it is.  But its effect is discontinuous.  Whether a  “Monday at 12:00” impuluse lasts for 7.494 days or 13.7781 days makes no difference: either way, it occurs on the first Monday and the second Monday, but not on the third Monday.

[Andrew Gelman]

Yes, but that's where my solution #1 or solution #2 could work.

[Phil Price]

In my implementation, which I have posted, duration is a continuous parameter (so is start time and end time).   This is your #1. This is already the way I do it. I'm not sure how many different ways I can say it.

Let me try it this way: Your solution #1 does not work.

As for #2...as I understand it, in order for this to work we need a function so steep that the smallest "continuous" change, at computer's precision, corresponds to an integer change in the parameter. Right?

That does seem like it might give the behavior I want, if it works at all. But then why not make this part of the language? It seems just incredibly kluldgy and fragile and machine-dependent and lots of other bad things.

I guess I should try it. I have to say, I'm pretty demoralized. It might be possible to trick Stan into fitting this model, but a lot of the appeal was the ability to write a model that expresses what I think is actually happening, and then fit it. If, instead, I'm going to have to find some delicately balanced model to trick Stan into fitting a type of model it's designed not to be able to fit, I sort of don't see the point.

[Bob Carpenter]

> On Dec 9, 2016, at 7:56 PM, Phil Price <pnp...@creekcats.com> wrote:
>

...

> I guess I should try it. I have to say, I'm pretty demoralized. It might be possible to trick Stan into fitting this model, but a lot of the appeal was the ability to write a model that expresses what I think is actually happening, and then fit it. If, instead, I'm going to have to find some delicately balanced model to trick Stan into fitting a type of model it's designed not to be able to fit, I sort of don't see the point.

Stan is hard to trick --- most of the attempts to do discrete parameters
other than marginalization will either just not work or they'll be very
slow to mix (the line between these isn't sharp).

We tried to build it to fit the largest class of continuous parameter
models that could be reliably fit. I don't think that includes the
models you're talking about because the discrete component is combinatorially
prohibitive to marginalize.

We might extend that goal to discrete parameters that can't be marginalized,
but that may just be a fool's errand as many of these systems are bound
to be truly intractable (in the sense of no polynomial algorithms) or
at best NP-complete (tractable unbounded parallel algorithms but no known
untractable algorithms and a very very strong suspicion none exist).

Let us know if you find something to fit this model as you conceive it.
You might want to look at PyMC3, which lets you write down a model with
mixed discrete and continuous parameters.

- Bob

[Andrew Gelman]

P

Yes, I understand your general concern.  I think the way Bob or Michael might put it is that purely discrete problems are essentially uncomputable because of the combinatorial explosion. So this leaves you with a few options:

1.  Set up the problem exactly as you did, accept the combinatorial explosion, and aim to get a reasonable solution, or some set of reasonable solutions.  This is an approximate-Bayes problem, and it represents an important general class of approximate-Bayes problems, and it's the sort of thing that Stan _should_ be able to do, but Stan can't yet do it.  So far with Stan we've focused on full Bayes, with approx algorithms as an afterthought.  But the machine learning world is the opposite:  they're mainly focused on approximate algorithms.  Those sorts of programs can require a bit of tuning, but this can be worth it if you can get good solutions.  Currently the machine learning programs seem to be focused on specific models.  It would be an interesting thing for someone to build out Stan to be able to attack this sort of problem.

2.  Back to my solution #1, which does not work.  Does it fail even if you simulate fake data with continuous parameters?  (As I wrote, "For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.")  I guess I'm still missing something here.  Yes, I understand that you are coding all these as continuous, but it also sounds like you're creating your data in such a way that the underlying parameters are discrete.  If you could have things start at 8:14pm on Tuesday, then wouldn't you be able to recover such things?

3.  In my solution #2, I'm assuming a function that's not so steep, that would allow the model to be fit but at the price of a bit of blurring of the inferences.  Such blurring is undesirable but it was my impression that you were already having a bit of blurring, a bit of difficulty identifying start and stop in the presence of noise.  So I was thinking that a little more blurring might not be so bad.

4.  As Bob wrote in his reply, you can't really trick Stan into fitting a discrete model that it can't fit.  My suggestions were more in the direction of making the problem easier, basically by blurring out the discrete modes.  I agree that, to the extent the real world is discrete, this is not an optimal approach, but it can work.  Imagine a spherical cow etc.

5.  Yet another approach is to keep the model discrete but to recognize that it has structure.  It's not a purely discrete problem because the weeks are ordered.  This would point toward some sort of custom solution to the problem, developing an algorithm that finds a good solution.  But that's a pain in the ass.  Pre-Stan, we wrote custom solutions all the time but now we prefer not to.

6.  Finally, maybe you can have it all (or most of it all) but first blurring the problem (following my ideas 1 and 2, which are items 2 and 3 in the above list) and then, if you want to work toward a discrete solution, using random draws from the continous-model posterior as starting points.

A

On Dec 9, 2016, at 7:56 PM, Phil Price <pnp...@creekcats.com> wrote:

In my implementation, which I have posted, duration is a continuous parameter (so is start time and end time).   This is your #1. This is already the way I do it. I'm not sure how many different ways I can say it.

Let me try it this way: Your solution #1 does not work.

As for #2...as I understand it, in order for this to work we need a function so steep that the smallest "continuous" change, at computer's precision, corresponds to an integer change in the parameter. Right?

That does seem like it might give the behavior I want, if it works at all. But then why not make this part of the language? It seems just incredibly kluldgy and fragile and machine-dependent and lots of other bad things.

I guess I should try it. I have to say, I'm pretty demoralized. It might be possible to trick Stan into fitting this model, but a lot of the appeal was the ability to write a model that expresses what I think is actually happening, and then fit it. If, instead, I'm going to have to find some delicately balanced model to trick Stan into fitting a type of model it's designed not to be able to fit, I sort of don't see the point.

Hmm, there are a couple of ways I can thinking of hacking this (or maybe "solving" rather than "hacking" if we can consider these changes as representing a true underlying model):

1.  Continuous-time parameter, continuous-time data:  Allow these weekly effects to start and end midweek.  For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.  With real data, if things actually only start at the beginning of each week, you might still be able to fit this model because the slop and uncertainty would allow a continuous fit.  Sure, it wouldn't be quite correct, but you already have enough noise that you're not always identifying the start and end _weeks_, so would it really be so horrible to allow continuous times?

2.  Continuous-time parameter, discrete-time data:  Again allow a continuous time for the start and end parameters.  But now match them to discrete weekly start and end via an S-curve.  For example, let's label the discrete points as integers.  Then if the continuous start parameter is x, the discrete starting point could be the integer closest to x + epsilon, where epsilon is a distribution such as a logistic or normal with center 0 and scale sigma, and sigma has some value such as 0.5 or 1:  this corresponds to an ordered logit or probit.  The idea is that sigma will be small enough that the start time parameter roughly localizes the actual start time, but large enough that the transition from one integer to the next is smooth, so there's no near-discontinuity that would jar the HMC.  You could implement this model with latent epsilon but I think it would probably work better to use the normal and logistic cdfs.  (In theory these go off to infinity but I'd just truncate to the nearest 3 or 4 discrete possibilities for the likelihood.)  Once you've fit the model, you can then, in the generated quantities block, use simple Bayes calculation to get inference for the epsilons and thus identify the actual starting and end points more precisely.

[Daniel Lee]

Just a quick thought regarding approximate algorithms: we can do some similar stuff in Stan by making approximate models. It's not quite the same, but you can write approximate models with single modes where advi or lbfgs will return something coherent. (so will mcmc, but on the approximate model)

Daniel

--

[Phillip Price]

On Dec 9, 2016, at 7:08 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

P

Yes, I understand your general concern.  I think the way Bob or Michael might put it is that purely discrete problems are essentially uncomputable because of the combinatorial explosion. So this leaves you with a few options:

1.  Set up the problem exactly as you did, accept the combinatorial explosion, and aim to get a reasonable solution, or some set of reasonable solutions.  This is an approximate-Bayes problem, and it represents an important general class of approximate-Bayes problems, and it's the sort of thing that Stan _should_ be able to do, but Stan can't yet do it.  So far with Stan we've focused on full Bayes, with approx algorithms as an afterthought.  But the machine learning world is the opposite:  they're mainly focused on approximate algorithms.  Those sorts of programs can require a bit of tuning, but this can be worth it if you can get good solutions.  Currently the machine learning programs seem to be focused on specific models.  It would be an interesting thing for someone to build out Stan to be able to attack this sort of problem.

I understand about the combinatorial explosion, but this is still a slightly unsatisfying situation. When a human looks at a time series that, starting on the 9th week, shows a 2-hour spike every Tuesday afternoon for the next 8 weeks, it’s easy to see what is happening. A computer can analyze chess moves a zillion times faster than me, so how come it’s paralyzed by the combinatorial explosion in this case and I’m not?  The answer, of course, is I’m not laboriously working through every single possibility. But why not?

Or think about some sort of modified traveling salesman problem, where I have some sort of noisy information on the total distance the salesman drove — maybe I know how much gas he used in his car. Let’s say he visits the 20 biggest cities in the U.S. and my estimate is that he drove 10% farther than the shortest possible route, plus or minus 8%, but no shorter than the shortest possible route of course.  99.99999% of possible routes are going to be a lot longer than that and will have probability very close to zero.  Is exhaustively exploring all of the routes really the best we can do? Clearly not. 

I understand Stan is built for a certain kind of model, and is just not the right tools for these kinds of jobs. But, as with your “R sucks” blog post, the fact that Stan is great at what it does seems to amplify my frustration when it won’t do this other thing I want it to do. 

In this case, in the continuous world Stan doesn’t force me to exhaustively sample all of continuous parameter space (something that can’t be done anyway, obviously!); by design, it spends most of its time in the high-probability regions of parameter space, and nearly none of its time where the probability is very very low. I would like to have the same capability in discrete space: hop around between discrete states, and spend most of its time where the probability is high. Indeed, if one of the tricks you suggested actually works (which I doubt) you could do this already in Stan: create a function that maps tiny “continuous” steps onto discrete state changes, and let Stan do its thing.

2.  Back to my solution #1, which does not work.  Does it fail even if you simulate fake data with continuous parameters?  (As I wrote, "For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.")  I guess I'm still missing something here.  Yes, I understand that you are coding all these as continuous, but it also sounds like you're creating your data in such a way that the underlying parameters are discrete.  If you could have things start at 8:14pm on Tuesday, then wouldn't you be able to recover such things?

Good question…I was thinking about this yesterday while exercising — it really is true that mulling things over while “not working” can be the best way to work —  and I think this approach will work, just not as I implemented it.  As with some of those morphing issues we never quite worked out, there might be a way to make this work as a fully continuous problem. After all, the actual electric load in the building really is a continuous function of time, even though it only gets measured every 15 minutes. 

This is going to be somewhat complicated, though, or at least I think so….Complicated to conceptualize and program, not necessarily to fit.  The issue is that there really are these discrete events, some of which recur daily and some weekly, and some only for limited duration, and they can overlap.  But yeah, I can  have weekly-recurring event 1 start at some time t1_0 and end at t1_1, and so on for other events, and since those times won’t match up to the meter times I’ll just do some interpolation at the ends… it all sounds workable, although I will have to be careful in order to define a fully specified model. 

For now I will do a toy problem and see if I can get it to work. I’m optimistic. Thanks for your help.

P

[Michael Betancourt]

In this case, in the continuous world Stan doesn’t force me to exhaustively sample all of continuous parameter space (something that can’t be done anyway, obviously!); by design, it spends most of its time in the high-probability regions of parameter space, and nearly none of its time where the probability is very very low. I would like to have the same capability in discrete space: hop around between discrete states, and spend most of its time where the probability is high

So would a huge number of scientists, statisticians, and

mathematicians.  We just don’t have any idea how to do

that for discrete spaces yet.  We’ve only recently figured

out how to do that for continuous spaces using gradient

information, but the correct generalization of “gradient” 

information to discrete spaces is not know.  This is why

we try to figure out ways of marginalizing down to 

continuous spaces and using the state of the art methods.

[Andrew Gelman]

On Dec 10, 2016, at 3:06 PM, Phillip Price <pnp...@creekcats.com> wrote:

On Dec 9, 2016, at 7:08 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

P

Yes, I understand your general concern.  I think the way Bob or Michael might put it is that purely discrete problems are essentially uncomputable because of the combinatorial explosion. So this leaves you with a few options:

1.  Set up the problem exactly as you did, accept the combinatorial explosion, and aim to get a reasonable solution, or some set of reasonable solutions.  This is an approximate-Bayes problem, and it represents an important general class of approximate-Bayes problems, and it's the sort of thing that Stan _should_ be able to do, but Stan can't yet do it.  So far with Stan we've focused on full Bayes, with approx algorithms as an afterthought.  But the machine learning world is the opposite:  they're mainly focused on approximate algorithms.  Those sorts of programs can require a bit of tuning, but this can be worth it if you can get good solutions.  Currently the machine learning programs seem to be focused on specific models.  It would be an interesting thing for someone to build out Stan to be able to attack this sort of problem.

I understand about the combinatorial explosion, but this is still a slightly unsatisfying situation. When a human looks at a time series that, starting on the 9th week, shows a 2-hour spike every Tuesday afternoon for the next 8 weeks, it’s easy to see what is happening. A computer can analyze chess moves a zillion times faster than me, so how come it’s paralyzed by the combinatorial explosion in this case and I’m not?  The answer, of course, is I’m not laboriously working through every single possibility. But why not?

Yes, this relates to my point 5 below, that although this problem is discrete, it's not quite a needle-in-a-haystack problem either.  It has structure.  We have good ways of thinking about structure with continuous models (normal distributions, hierarchical models, etc), but not so much with discrete models.    If all discrete possibilities were truly possible, then there would be a combinatorial explosion and this problem would not be possible to solve "by eye," either.  But we actually have lots of prior information, which is not so easy for us to include in the model. 

Or think about some sort of modified traveling salesman problem, where I have some sort of noisy information on the total distance the salesman drove — maybe I know how much gas he used in his car. Let’s say he visits the 20 biggest cities in the U.S. and my estimate is that he drove 10% farther than the shortest possible route, plus or minus 8%, but no shorter than the shortest possible route of course.  99.99999% of possible routes are going to be a lot longer than that and will have probability very close to zero.  Is exhaustively exploring all of the routes really the best we can do? Clearly not. 

I agree.  That's what I was talking about in my point 1 above, that there should be some general ways of getting good approximate solutions, stepping back from our traditional goal of full Bayesian inference.  I'd love it if we could do this in Stan, and I think it's possible.  It would take some work, thought!

I understand Stan is built for a certain kind of model, and is just not the right tools for these kinds of jobs. But, as with your “R sucks” blog post, the fact that Stan is great at what it does seems to amplify my frustration when it won’t do this other thing I want it to do. 

In this case, in the continuous world Stan doesn’t force me to exhaustively sample all of continuous parameter space (something that can’t be done anyway, obviously!); by design, it spends most of its time in the high-probability regions of parameter space, and nearly none of its time where the probability is very very low. I would like to have the same capability in discrete space: hop around between discrete states, and spend most of its time where the probability is high.

Yes, I agree that it would be good to be able to do this.

2.  Back to my solution #1, which does not work.  Does it fail even if you simulate fake data with continuous parameters?  (As I wrote, "For example, suppose that the "start time" is 2.39, which corresponds, say, to 8:14pm on Tuesday.  Then the change would start at that time.  With fake data, you could actually do this by picking the start time from some distribution and then starting things at whatever continuous time comes out.")  I guess I'm still missing something here.  Yes, I understand that you are coding all these as continuous, but it also sounds like you're creating your data in such a way that the underlying parameters are discrete.  If you could have things start at 8:14pm on Tuesday, then wouldn't you be able to recover such things?

Good question…I was thinking about this yesterday while exercising — it really is true that mulling things over while “not working” can be the best way to work —  and I think this approach will work, just not as I implemented it.  As with some of those morphing issues we never quite worked out, there might be a way to make this work as a fully continuous problem. After all, the actual electric load in the building really is a continuous function of time, even though it only gets measured every 15 minutes. 

This is going to be somewhat complicated, though, or at least I think so….Complicated to conceptualize and program, not necessarily to fit.  The issue is that there really are these discrete events, some of which recur daily and some weekly, and some only for limited duration, and they can overlap.  But yeah, I can  have weekly-recurring event 1 start at some time t1_0 and end at t1_1, and so on for other events, and since those times won’t match up to the meter times I’ll just do some interpolation at the ends… it all sounds workable, although I will have to be careful in order to define a fully specified model. 

3.  In my solution #2, I'm assuming a function that's not so steep, that would allow the model to be fit but at the price of a bit of blurring of the inferences.  Such blurring is undesirable but it was my impression that you were already having a bit of blurring, a bit of difficulty identifying start and stop in the presence of noise.  So I was thinking that a little more blurring might not be so bad.

4.  As Bob wrote in his reply, you can't really trick Stan into fitting a discrete model that it can't fit.  My suggestions were more in the direction of making the problem easier, basically by blurring out the discrete modes.  I agree that, to the extent the real world is discrete, this is not an optimal approach, but it can work.  Imagine a spherical cow etc.

5.  Yet another approach is to keep the model discrete but to recognize that it has structure.  It's not a purely discrete problem because the weeks are ordered.  This would point toward some sort of custom solution to the problem, developing an algorithm that finds a good solution.  But that's a pain in the ass.  Pre-Stan, we wrote custom solutions all the time but now we prefer not to.

6.  Finally, maybe you can have it all (or most of it all) but first blurring the problem (following my ideas 1 and 2, which are items 2 and 3 in the above list) and then, if you want to work toward a discrete solution, using random draws from the continous-model posterior as starting points.

[Phillip Price]

The good news is, I have found a solution to my problem: Now a key parameter (“week number”) is treated as a truly continuous variable rather than a discrete one. 

The bad news is, the program is still really, really slow (as in: three hours to reach reasonable convergence on a fairly simple test problem). I don’t understand why. 

Here’s how I made the problem continuous (skip the stuff between the === below if you don’t need the details)

===

I have a coefficient betaTOD[ week, TOD] that quantifies the “time of day (TOD) effect” in a given week. If, in week 4, at 9:00 every day the electric load is 2 kW, then betaTOD[ 4, 9] = 2

I write this coefficient as a sum: betaTOD[ week, TOD] =  betaTODinit [ TOD ] + phiTOD[ week, TOD] 

betaTODinit is sampled from a prior

phiTOD is a step function that starts in one week and ends in another. It has a start week defined by weekStartTOD, a magnitude defined by deltaTOD, and a duration defined by tauTOD. 

I want to be able to have “impulses” that last for some number of weeks at a given time of day. At 9 am, if the load is 1 kW in weeks 1 and 2 and 3, and then goes up to 2 kW for weeks 4 and 5 and 6, and then drops back to 1 kW, I want to capture this as an impulse that has a height of 1 kW, a duration of 3 weeks, and starts in week 4. 

The problem is that, as described, this is a discrete problem even if I make the start time and the duration continuous: whether the “9 am” effect starts when week = 3.001 or 3.999 makes no difference: either way, the impulse doesn’t apply in week 3 but it does in week 4.  (In the full model I want to fit, there’s a reason I don’t start an impulse mid-week; otherwise, it would make sense for week = 3.5 to mean the impulse starts on Wednesday. But for this toy problem, an impulse either applies to the whole week at a given time of day, or it doesn’t apply at all that week). 

===

My solution, which I was SURE was going to work, was to change the definition of the model so that a non-integer start time has an effect on the predictions. Suppose the candidate start time is 3.001. Previously, what this meant is that the first week in which the effect occurs is week 4, which the first week after the start of the impulse.  What I did was change this so that now some fraction of the effect shows up _before_ the start time. What fraction? 

f =  (integer start week - continuous start week) So, if the candidate start time is week = 3.001, this generates an integer start week of 4, and we get f= 0.999. So even though the “start time” is just AFTER 3.0, we get 99.9% of the effect during week 3. In contrast, if week = 3.999, we only get 0.1% of the effect in the week 3. 

The effect of all of this is that the program can converge on the right answer along a continuous path.  If it tries a start time of 3.001 for an impulse that actually started in week 4, it will apply a lot of that impulse to week 3 and thus overestimate that week (assuming the impulse is positive). It can gradually work its way towards a start time of 4.0 by going downhill (towards lower error) all the way. Of course Stan isn’t trying to just run downhill, it’s sampling, but the point is it now has a nicely behaved continuous model for doing this. And it works!

But it works…so…very….slowly. 

I’m doing an 11-week test problem, 24 times of day, 5 days per week. I have only one “impulse” for each time of day (e.g. there can be one start time and duration and magnitude for 1:00, and a different start time and duration and magnitude at 2:00, and so on).  So, 24 impulses, each defined by 3 adjustable parameters, plus the 24 parameters that define the initial state…altogether that’s about 100 parameters. How can it take hours?

I’ve taken a look at the “pairs” plot and it looks like the model is spending all of its time in a few regions of parameter space that aren’t near the maximum. See the attached plot. Here I did 800 iterations and as you can see t

. The “right answer” is:

betaTODinit[10] = 1.0,  (This is the magnitude of the “hour 10” effect in week 1)

deltaTOD[1,10] = 1.0,  (This is the magnitude of the impulse)

weekStartTOD[1,10] = 3.0, (This is the start time of the impulse)

sigma[10] = 0.05 (This is the amount of noise added to the signal, which would otherwise be a square wave)

If it finds the right parameters for these and for all of the other parameters in the model — everywhere there’s a 10 in the list above, there can be any number from 1 to 24 — it can have its predictions match the data almost exactly, the only difference is iid normal noise with sd = 0.05.  But unless I cook the books by giving it very tight priors around the right answer, it doesn’t find it. 

Basically it looks like its step size is way too small. But I’ve tried increasing adapt_delta to 0.95 (from the default of 0.8) and this doesn’t seem to help. 

If I print out some values for debugging, it typically makes really really tiny steps (near some incorrect answer) until suddenly it takes a big jump towards a somewhat better answer and then takes really really tiny steps again. You can see this in the pairs() plot. As I’m writing this, I’m running it while I watch it print out the maximum residual. The model is trying to pull this towards zero: y ~ N( yhat, sigma) and sigma ~ half-normal(0, sigmaMax).  Anyway, it spend hundreds of iterations with a big max residual, around 2.2, changing only the fourth or fifth digit. Then suddenly it ran down to about 1.3 and then stalled there. Then down to 0.9 and stalled there….and that chain just ended at 1000 iterations and now it is doing another one, starting with the max residual around 1.45 again only changing the fourth digit as I write this…whoah, except it just jumped down to 0.76…and now back up to 0.93, where it is changing only in the fourth decimal place. 

You’ve said in the past that slow convergence usually means a bug in the model, and of course that’s possible here…but the fact that I can give it tight priors around the right answer and it agrees that yep, that’s the right answer…that makes me think there’s not a bug.

There is inevitably some correlation between parameter estimates — for example,  if it has found the “right” end time but its start time is too early, it needs to both move the start time later and decrease the duration in order to make things better — but even if it were just picking completely at random, half the time it tries a later start time it should try a lower duration. Surely it should not take thousands and thousands of iterations to make progress.

I’m baffled.

Ideas welcome!

P

[Andrew Gelman]

Hmmm, hard to do this over email . . . my quick thought is the generic recommendation to strip down the model and see if everything's ok on a simpler version. From the pairs plot it looks like things are multimodal which makes me really tempted to use this as an example for my new adaptive simulated tempering algorithm. But I wonder whether there's a way to make the model even more continuous in some way. Also I wonder if there's some real nonidentifiability because of the starting and stopping model.

Just for example--and plz recognize that I'm making it up becuase I'm not trying to follow all your details--suppose you have something that starts at time 1.3 and ends at time 10.2. Then couldn't it also be something that starts at 1.3 and ends at X, plus another something that starts at X and ends at 10.2? You get the idea. If there are lots of these impulses and they can adjoin or overlap, then there could be lots of stories that will legitimately fit the data. Which suggests maybe some stronger priors are needed about the legnths of these impulses, or the time between them, or something like that.

A

> <pairsPlot.png>

[Phillip Price]

> On Dec 15, 2016, at 10:48 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:
>
> Hmmm, hard to do this over email . . . my quick thought is the generic recommendation to strip down the model and see if everything's ok on a simpler version. From the pairs plot it looks like things are multimodal which makes me really tempted to use this as an example for my new adaptive simulated tempering algorithm. But I wonder whether there's a way to make the model even more continuous in some way. Also I wonder if there's some real nonidentifiability because of the starting and stopping model.
>
> Just for example--and plz recognize that I'm making it up becuase I'm not trying to follow all your details--suppose you have something that starts at time 1.3 and ends at time 10.2. Then couldn't it also be something that starts at 1.3 and ends at X, plus another something that starts at X and ends at 10.2? You get the idea. If there are lots of these impulses and they can adjoin or overlap, then there could be lots of stories that will legitimately fit the data. Which suggests maybe some stronger priors are needed about the legnths of these impulses, or the time between them, or something like that.
>

Right now I am working with a single impulse in both the data and the model, so this is not what is happening.

And actually the situation is worse than I thought. I will boil it down to an even simpler test model.

P