Optimisation issues with two time series involving transfer function

[Tom Elliott]

Hello everyone,

I'm working on my honours project at the University of Auckland with Dr. Brendon Brewer. It involves modelling reverberation mapping data for Active Galactic Nuclei. More simply, the data are two time series. We need to first fit an AR(1) model to the first time series (Y1) that only has partial observations (with known measurement error), and then use this to model the second time series (Y2) which follows the first through a transfer function, whose parameter of interest is tau, the lag between the first and second time series.

The formula for the first time series is:
Y1[t] = mu + alpha * (Y1[t-1] - mu), a typical AR(1) model, with the observations y1(t) ~ Normal(Y1[t], beta^2). For the data we are using approximately 1000 points in the full Y1 time series, and we have roughly 10% of them observed (the reason for needing so many points is that we are approximating a continuous autoregressive CAR(1) model). If we apply this model to the first time series alone, STAN converges nicely and estimates the parameters very well.

The second series is where we start running into issues. We are assuming a uniform transfer function between two values, a and b, so that the value of the second series at time t, Y2[t], is the mean of the first time series between times t-a and t-b, with some other constants for variability and shift (eta and C).
The equation is:
Y2[t] = eta * ( (1 / N) * SUM(Y[t - a] : Y[t - b]) + C ), where N is the number of values in the summed interval.

The problem we are having is summing over the correct range of Y1. STAN doesn't allow subsets (as far as I can tell), so I have used an if() ... control flow to go through every Y1 value and check if the time value, t, is in the required interval.

This method is working, however the model runs quite slowly, and takes a couple of hours to analyse the data from one object, and in the end we will need to apply the model to many objects (~100).

I have attached a working STAN model, and would like to know if you have any ideas or techniques that will help achieve convergence faster. The autocorrelation is very high, however because the chains take so long to run I am hesitant to simply use a large trimming interval, so any ways to try and minimise this would very much be appreciated.

I have also included an R script (example.R) that includes the type of data we will be using, incase you want to test any ideas you have.

Thanks in advance for any help you can offer,

Tom Elliott.

[Ben Goodrich]

On Friday, July 19, 2013 12:06:53 AM UTC-4, Tom Elliott wrote:

The second series is where we start running into issues. We are assuming a uniform transfer function between two values, a and b, so that the value of the second series at time t, Y2[t], is the mean of the first time series between times t-a and t-b, with some other constants for variability and shift (eta and C).
The equation is:
Y2[t] = eta * ( (1 / N) * SUM(Y[t - a] : Y[t - b]) + C ), where N is the number of values in the summed interval.

The problem we are having is summing over the correct range of Y1. STAN doesn't allow subsets (as far as I can tell), so I have used an if() ... control flow to go through every Y1 value and check if the time value, t, is in the required interval.

The next version of Stan has a segment() function that will make your life easier. The way you have done it now is fine and doesn't really slow the execution time down much. More likely, the more complicated model makes it harder for the MCMC algorithm to adapt.

I am running your example.R now, but I guess it will take me overnight. If you have your stan.fit object saved, I would be interested to see the plot produced by

pairs(stan.fit, pars = c("alpha", "beta", "mu", "b", "k", "eta", "C", "lp__"))

which should illuminate pairs of parameters that are highly correlated.

Ben

[Tom Elliott]

Hi Ben, here's the ouput from the pairs plot.

[Ben Goodrich]

On Friday, July 19, 2013 1:08:02 AM UTC-4, Tom Elliott wrote:

Hi Ben, here's the ouput from the pairs plot.

This is a great example of why I like the pairs() plot. So, the issues that I can see are

1. The distribution of alpha is pushed up against its upper bound of 1
2. mu is not mixing very well
3. The distribution of b is possibly pushed up against its lower bound of 0
4. Something weird is happening with eta and C
5. The row / column with lp__ in it does not look like it should under a well-performing model (roughly Gaussian with density concentrated near the mode)

You might compare with the corresponding plot from the simpler model you said was working well. Unfortunately, I am about to fall asleep and can't suggest any solutions right now. Maybe Bob can pick it up for the next few hours.

Ben

[Michael Betancourt]

There's nothing weird here -- with Y2[t] well-constrained we would expect eta \approx 1 / C which is exactly what is seen in the plot. My guess is that the prior on C is too weak and it's growing to the point of dominating the SUM(Y[t - a]: Y[t - b]) / N contribution. A stronger prior pushing C towards zero may help.

[Michael Betancourt]

Actually, it should be easier to just factor C out, i.e. Y2[t] = eta * ( (1 / N) * SUM(Y[t - a] : Y[t - b])) + C;

[Ben Goodrich]

On Friday, July 19, 2013 4:39:42 AM UTC-4, Michael Betancourt wrote:

There's nothing weird here -- with Y2[t] well-constrained we would expect eta \approx 1 / C which is exactly what is seen in the plot.  My guess is that the prior on C is too weak and it's growing to the point of dominating the SUM(Y[t - a]: Y[t - b]) / N contribution.  A stronger prior pushing C towards zero may help.

Ah, now that I am fully awake I can see that makes sense. But if so, then a model {} block that has independent diffuse priors on C and eta doesn't seem quite right, so I would agree that a different prior on C may help.

The .stan file also comments that a

alpha ~ beta(20,1)

prior is necessary to prevent the chains from getting stuck at unlikely values. Can you say more about this? A prior as strong as that may be appropriate in your field, but it is going to make it difficult to estimate the other parameters precisely.

Ben

[Tom Elliott]

The reason for such a strong prior on alpha is that when we used a uniform, it tended to favor small alpha values and essentially fit a straight line through the points and and let beta be very big, which is complete nonsense.

I tried giving C a much narrower prior and this seems to have stopped it from roaming around so much and instead hang around 5 which is what it is when we use the full arp151 data set. Other than this I'm not sure on an alternative prior for C.

[Michael Betancourt]

I real issue isn't C but rather the product eta * C, which isn't fully specified no matter how tight your data. Stronger priors can help, as you've seen, as would removing the product by moving the offset to outside the scaling.

On Jul 19, 2013, at 10:11 PM, Tom Elliott <tomell...@gmail.com> wrote:

> The reason for such a strong prior on alpha is that when we used a uniform, it tended to favor small alpha values and essentially fit a straight line through the points and and let beta be very big, which is complete nonsense.
>
> I tried giving C a much narrower prior and this seems to have stopped it from roaming around so much and instead hang around 5 which is what it is when we use the full arp151 data set. Other than this I'm not sure on an alternative prior for C.
>

> --
> 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/groups/opt_out.
>
>