Trouble estimating exponential decay rate from time series

Modeling
,
1

Hey all!

I’m trying to estimate the exponential rate of decay from a time series, and am running into some problems.

I’ve been using Stan for a few years now, mostly to fit ODEs, so I’m sure I’m doing something incredibly obviously wrong, but would appreciate help nonetheless!

My data look like this:

y_samp_repro <- c(49807, 31437, 17739, 10739, 6499, 3836, 2201,  1240, 823, 436, 259, 179, 91, 61, 46, 17, 16, 7, 4, 3, 2, 1)

I’m using this paper as a guide, so started out writing the Stan model as follows:

data {
  int LENGTH;
  vector[LENGTH] Y;
}
parameters {
  real lambda;
}
model {
  real alpha;
  real beta;
  alpha = 1;
  beta = 1;
  lambda ~ gamma(alpha, beta);
  Y ~ exponential(lambda);
}
generated quantities {
  real pred;
  pred = exponential_rng(lambda);
}

In the paper, they’re looking to estimate the rate of exponential growth, but I’m trying to estimate the rate of decay, which might be where I’m running into issues?

When I fit the model:

model <- stan_model("exp-decay.stan")
fit <- sampling(object = model,
                data = list(Y = y_samp_repro, LENGTH = length(y_samp_repro)),
                warmup = 500,
                iter = 1000,
                chains = 2,
                cores = 2,
                seed = 123); print(fit)

I get an output, but it’s not very reasonable:

Inference for Stan model: anon_model.
2 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=1000.

          mean se_mean     sd    2.5%     25%     50%     75%    97.5% n_eff Rhat
lambda    0.00    0.00    0.0    0.00    0.00    0.00    0.00     0.00   400    1
pred   5867.79  191.86 6059.9  167.91 1765.62 4029.63 7692.35 20953.34   998    1
lp__   -212.78    0.04    0.7 -214.89 -212.99 -212.52 -212.32  -212.27   327    1

Samples were drawn using NUTS(diag_e) at Mon Jul 17 17:34:11 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Again, I’m sure the problem is really, really obvious (and suspect it has to do with looking at decay rather than growth), but I’m just not really sure where to start, so would appreciate any help at all!

Thank you so much for any help!

2

Hi @MadelineJC , I might be misunderstanding something here, but your data are bivariate (you write that this is a time series, so I’m thinking that y_samp_repro is some response observed across time), right? I plotted your data on semilog axes, and they look linear, presuming that the time is evenly spaced. This wouldn’t fit with the exponential distribution, which is univariate.

Exponential decline in time series data leads me to think you would want a model of the form y=e^{-k\cdot t}.

2 Likes
3

I agree that it is most likely the case here.
So maybe given how the data looks it could be something like: Poisson(y | e^{−kt}).

1 Like
4

Hi AWoodward! I’m sure the misunderstanding is mine… y_samp_repro is a response observed overtime, generated from an exponential decay function (as written in your response). Yes, when you plot the logged data, they are linear, but I assumed that because the data were generated using an exponential decay function, I should be trying to estimate the rate of decay (k) by fitting an exponential model. I’ll try playing around with the Poisson distribution, as suggested below! Thank you for your response!

5

Thanks for your response! I’ll trying playing around with the Poisson distribution!

6

Take it with a grain of salt though. I just assumed by the data points that you might deal with counts. It could still be i. e. radioactive counts but you should check what fits best ;)