# Marginalising 2 parameters: Poisson change point model

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

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.

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

Could you post the working code? Many thanks.