Change Point Model

[Antonio Martini]

I was reading chapter '14.2 Change Point Models' of the Stan manual, specifically the example and model related to the coal mining disasters.

Could the conditional operator be replaced by a continuous step mixing function in order to simplify the model?

I mean  replacing:

 Dt ~ Poisson(t < s ? e : l)  

With something like:

  alpha = 1.0 / (1 + exp(s - t - 0.5)^15);                     // step function modeling rate change

  Dt ~ Poisson( e*(1-alpha) + l*alpha );

We can see that at the change point alpha is :

1.0 / (1 + exp(0 - 0.5)^15) = 0.9994472 

and 

1.0 / (1 + exp(1 - 0.5)^15) =  0.0005527786

I have implemented the example in full and the results are consistent with all the other solutions I have seen. Full code below.

Is this a valid way of modelling this problem in Stan?

Thanks,

Antonio

------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# coal mining disasters change point test

require(Rcpp)

require(rstan)

stanmodel <- '

data {

  // signal data

  int<lower=1>          N;

  int                       y[N];          // value

}

parameters {

  real<lower=1>             event;       // change point day

  real<lower=0>             lambda[2]; // rates before and after change point

}

model {

  event  ~ uniform(1,112); // [1, max(day)]

  lambda ~ uniform(0,10);  // [1, max(y)]

  for (i in 1:N)

  {

    real alpha;

    real beta;

    alpha = 1.0 / (1 + exp(event - i - 0.5)^15);                     // step function modeling change point

    beta  = lambda[1] *(1-alpha) +  lambda[2]*alpha;          // rates mixing

    y[i] ~ poisson(beta);

  }

}

'

# data

y <- c(4,5,4,1,0,4,3,4,0,6,

        3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,1,5,5,3,4,2,5,2,2,3,4,

        2,1,3,2,1,1,1,1,1,3,0,0,1,0,1,1,0,0,3,1,0,3,2,2,0,1,

        1,1,0,1,0,1,0,0,0,2,1,0,0,0,1,1,0,2,2,3,1,1,2,1,1,1,

        1,2,4,2,0,0,0,1,4,0,0,0,1,0,0,0,0,0,1,0,0,1,0,0)

 

N = length(y)

mdl <- stan(model_code=stanmodel, data= c("N", "y" ), iter=4000, chains=8, cores=8, control = list(adapt_delta = 0.95) )

# extra results and plot

t <- extract(mdl)$event  # change point distribution'

day <- 1:N

par(mfrow=c(2,1))

hist(t, col='lightBlue', freq=F, xlab="year", ylab="P(t)", main="", breaks=40)

abline(v=mean(t), col='red')

plot(y, cex=0.8, xlab='year', ylab='d')

lines(y)

abline(v=mean(t), col='red')

[Bob Carpenter]

Yes, and in fact Andrew's often suggesting this kind of
underlying continuous model. The trick is making everything
continuously differentiable like you have.

You may be interested in checking out is Poisson
processes, which provide natural geo-spatial models for
count data.

Does anyone have any examples of those coded up in Stan
they can share?

- Bob

> On Feb 24, 2017, at 10:50 AM, Antonio Martini <antonio...@hotmail.com> wrote:
>
>
> I was reading chapter '14.2 Change Point Models' of the Stan manual, specifically the example and model related to the coal mining disasters.
>
> Could the conditional operator be replaced by a continuous step mixing function in order to simplify the model?
> I mean replacing:
>
> Dt ~ Poisson(t < s ? e : l)
>
> With something like:
>

> alpha = 1.0 / (1 + exp(change_point - t - 0.5)^15); // step function modeling change point

> Dt ~ Poisson( e*(1-alpha) + l*alpha );
>
> We can see that at the change point alpha is :
>

> 1.0 / (1 + exp(- 0.5)^15) = 0.9994472

> hist(t, col='lightBlue', freq=F, xlab="day", ylab="P(t)", main="", breaks=40)
> abline(v=mean(t), col='red')
> plot(y, cex=0.8, xlab='day', ylab='orders')

> lines(y)
> abline(v=mean(t), col='red')
>

> --
> 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.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Antonio Martini]

Great to hear that as it saves a lot of extra work. 

Thanks! Antonio