Marginalising 2 parameters: Poisson change point model

Modeling
1

Hello

I am trying to write a poisson arrival model that incorporates 2 rate changes. Unlike the change point model in the manual (8.2), my change points are in continuous time. For the single rate change model I managed to marginalise the parameter using integrals and achieved a similar result to that of the manual’s example. Now with the 2 rate change model I am getting the fatal error of ‘Log probability evaluates to log(0)’. I was wondering if anyone has experience with a similar model.

This is my simulated data I am using to test the model
simuldata <- rep(0,23)
simuldata[1:7] <- rexp(7,0.001)
simuldata[8:15] <- rexp(8,0.002)
simuldata[16:23] <- rexp(8,0.003)
D = cumsum(simuldata)
T = cumsum(simuldata)[length(simuldata)]+100

data {
  real<lower=0> T;                //total observation time
  int<lower=1> N;                 //number of observations
  vector[N] D;                    //arrival time data
}

transformed data{
  vector[N+2] DAT;                
  DAT[1] = 0;
  DAT[2:(N+1)] = D;
  DAT[N+2] = T;                   //appends bracket of total observation time around data
}

parameters {
  real<lower=0> lambda1;          //first rate
  real<lower=0> lambda2;          //second rate
  real<lower=0> lambda3;          //third rate
}
transformed parameters {
  matrix[(N+1),(N+1)] lp3;
  for (i in 1:(N+1)) {
    lp3[i,i] = log((lambda1/lambda3)^(i-1)*((exp(-(lambda1-lambda2)*DAT[i+1])-exp(-(lambda1-lambda2)*DAT[i]))*exp(-(lambda2-lambda3)*DAT[i+1])/(lambda2-lambda1) - 
                                      (exp(-(lambda1-lambda3)*DAT[i+1])-exp(-(lambda1-lambda3)*DAT[i]))/(lambda3-lambda1))/(lambda3-lambda2));
  }
  for (i in 1:N) {
    for (j in (i+1):(N+1)) {
      lp3[i,j]= log(((lambda1/lambda2)^(i-1))*(exp(-(lambda1-lambda2)*DAT[i+1])-exp(-(lambda1-lambda2)*DAT[i]))*
                   ((lambda2/lambda3)^(j-1))*(exp(-(lambda2-lambda3)*DAT[j+1])-exp(-(lambda2-lambda3)*DAT[j]))/((lambda2-lambda1)*(lambda3-lambda2)));
    }
  }
}
model {
  target += N*log(lambda3) - lambda3*T - log((T^2)/2);
  target += log_sum_exp(to_vector(lp3));
  lambda1 ~ exponential(580);
  lambda2 ~ exponential(580);
  lambda3 ~ exponential(580);
}

If needed I will provide the derivation for the calculations, however I’m pretty certain they are correct as I can run the loops in R and the values are reasonable.

Thanks

1 Like
2

Here’s some generic advice to start with. Why don’t you write the complicated bits in the model block as functions in the functions block? This would allow you to export the exact same functions you’re using in your Stan program to R (via expose_stan_functions) and then test in R. I’ve found that often I think I’m doing the same thing in R and Stan but I’m not.

3

Thanks for the advice, I think if I had done this originally I would have avoided my problems. I was ignoring the unassigned values of the matrix, which meant they were included in the log_sum_exp at the end, hence the log(0) error

4

Could you post the working code? Many thanks.