Using a continuous distribution to model a time series

[Jan]

Hi,

I was wondering if one could use the distributions included in Stan to model a time series. E.g. a process that starts at zero then rises quickly and returns to baseline more slowly could be approximated by a gamma density (or mixture of gamma densities).

Similar to linear models (y_hat <- a + b*x) one would declare this relationship in the transformed parameters block, I assume. I couldn't find a good example, except for  example-models/bugs_examples/vol2/pines/pines-2.stan which uses the normal_log function in the transformed block.

As described in the attached files, I simulate a gamma density-based time series

x <- 1:10
y <- dgamma(x, shape=alpha, scale=1/beta)
y <- y + rnorm(length(x), 0, sigma_y)

and then try to recover the gamma parameters with stan

data {
  int<lower=0> N;  // nr samples
  real y[N];       // y values
  real x[N];       // x values
}
parameters {
  real<lower=0> alpha;
  real<lower=0> beta;
  real<lower=0,upper=100> sigma;
}
transformed parameters {
  real y_hat[N];
  for(n in 1:N)
    y_hat[n] <- gamma_log(x[n], alpha, beta);
}
model {
  for(n in 1:N)
    y[n] ~ normal(y_hat[n], sigma);
}

This sort-of works, but estimated parameters become non-sensical when adding more terms to the data and model, e.g.

y <- a + b * y + rnorm(length(x), 0, sigma_y)

I was wondering if I should abandon this approach in favour of a better way to parameterize smooth time series?

Cheers.

[Bob Carpenter]

You are simulating y on the linear scale and then
using the log in Stan. To match your simulation,
you need to do this:

y_hat[n] <- exp(gamma_log(x[n], alpha, beta));

I'll let others comment on the wisdom of this kind of model,
which is going to depend on how your data's generated.
And then you're going to fit

We also don't recommend the uniform(0, 100) priors --- see
the manual for current recommendations and Andrew's to-appear
Wiki. And then you want to vectorize this for efficiency:

> for(n in 1:N)
> y[n] ~ normal(y_hat[n], sigma);

- Bob

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

[Jan]

Thanks for the quick reply, Bob.

I tried adding the exp, but that causes a bad fit and a lot of errors:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Exception thrown at line 17:

Error in function boost::math::lgamma<double>(double): numeric overflow

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

I'm not sure how wrapping the gamma_log with exp can make the gamma_log overflow?

Thanks.

[Bob Carpenter]

Did you manage to get it to run until convergence? Did you
keep getting these errors? You can try lowering initial step size.

gamma_log() returns log(gamma(...)), so without the exp
to convert it back, you're not going to get sensible results.

You could move to modeling on the log scale.

- Bob

[Krzysztof Sakrejda]

Why not try some more flexible model?  Splines come to mind.  Krzysztof

[Jan]

The errors were cause by too large a values. Putting a more restrictive prior on the parameters of the gamma density helped, e.g.

model {
 y ~ normal(y_hat, sigma);
 alpha ~ cauchy(0,10);  // stay way below 100 here
 beta  ~ cauchy(0,10);  // dito
 sigma ~ cauchy(0,10);  // not necessary
}

The model now converges even with 2k iterations. So restricting it to positive uniform is not enough. Playing around with stan(..., control=list(stepsize=0.1)) did not have the same results as more restrictive priors.

I agree, Krzysztof, that a spline model would be more useful for a more variable time course sampled at many time points. However, here I'm looking at a special case where the response of a system is very simple and the less parameters I have to estimate from the noisy and coarsely sampled data the better. A sensible knot placement for a spine model in such a case might be difficult to achieve. But correct me if I'm wrong.

Cheers.

[Bob Carpenter]

Modifying step size just reduces numerical errors --- you
shouldn't get a different answer if everything converges without
exceptions during sampling.

We've tended to jump straight from parametric models to Gaussian
processes --- if someone has some nice spline models they want
to share, I'd love to see them (sorry if they've already been
posted and I missed them) --- then I might be able to finally
understand spline models, at which point I can put some examples
in the manual.

- Bob

> On Aug 25, 2015, at 9:50 PM, Jan <aaka...@googlemail.com> wrote:
>
> The errors were cause by too large a values. Putting a more restrictive prior on the parameters of the gamma density helped, e.g.
>
> model {
> y ~ normal(y_hat, sigma);
> alpha ~ cauchy(0,10); // stay way below 100 here
> beta ~ cauchy(0,10); // dito
> sigma ~ cauchy(0,10); // not necessary
> }
>
> The model now converges even with 2k iterations. So restricting it to positive uniform is not enough. Playing around with stan(..., control=list(stepsize=0.1)) did not have the same results as more restrictive priors.
>

> I agree, Krzysztof, that a spline model would be more useful fore a more variable time course sampled at many time points. However, here I'm looking at a special case where the response of a system is very simple and the less parameters I have to estimate from the noisy and coarsely sampled data the better. A sensible knot placement for a spine model in such a case might be difficult to achieve. But correct me if I'm wrong.
>
> Cheers.
>
>
>
>
>

[Jan]

I might end up looking into GP or splines after all. It turns out that extending the model with intercept and scaling factors a, b

y_hat[n] <- a + b * exp(gamma_log(x[n], alpha, beta));

introduces horrible correlations between the four parameters. No wonder I had problems with convergence. My gut feeling is that such a model would need to be re-paramterized in a clever way to avoid these problems. But I'll try normalizing the input data first.

[Krzysztof Sakrejda]

On Wednesday, August 26, 2015 at 12:07:01 AM UTC-4, Bob Carpenter wrote:

Modifying step size just reduces numerical errors --- you
shouldn't get a different answer if everything converges without
exceptions during sampling.

We've tended to jump straight from parametric models to Gaussian
processes --- if someone has some nice spline models they want
to share, I'd love to see them (sorry if they've already been
posted and I missed them) --- then I might be able to finally
understand spline models, at which point I can put some examples
in the manual.  

What kind of use cases would be relevant for the manual?  I was going to post something
here but all my data are a little exotic.

Krzysztof

 

[Bob Carpenter]

Just a simple example of how to fit a spline. As a computer
scientist, I prefer the simplest model possible to illustrate
a feature rather than something more realistic.

- Bob

[Michael Betancourt]

As in “fit 1D data” simple. Say simulated from a high-order polynomial.

[Bob Carpenter]

Sure. That'd be a fine example. I don't know how to do the
predictions with splines. I do think I finally understand
the basics of GPs.

- Bob

[Krzysztof Sakrejda]

Ok, I see what you're going for here.  I don't have anything like that on hand but I'll have something to contribute after the plumbing the sparse stuff through stan.  Krzysztof