#Bayes 1 model with up to 4 joinpoints
model{
  for(i in 1:n){
    Obs[i]~dpois(mean[i])

    #Modelling of the expected cases for every period
    #the maximum and minimum bounds are imposed to avoid numerical overflow.
    #x1 and x2 are the initial and final years of study.
    log(mean[i])<-log(Pob[i])+alpha+beta0*(anyo[i]-(x1+x2)/2)+
    min(max(delta[1]*beta[1]*BaseJoin[i,1]+
    delta[2]*beta[2]*BaseJoin[i,2]+
    delta[3]*beta[3]*BaseJoin[i,3]+
    delta[4]*beta[4]*BaseJoin[i,4],-10),10)
    
    #estimated rate at every period
    rate[i]<-exp(alpha+beta0*(anyo[i]-(x1+x2)/2)+min(max(delta[1]*beta[1]*BaseJoin[i,1]+
    delta[2]*beta[2]*BaseJoin[i,2]+delta[3]*beta[3]*BaseJoin[i,3]+delta[4]*beta[4]*BaseJoin[i,4],-10),10))
  }

  #APC at every period
  for(i in 1:(n-1)){APC[i]<-100*(rate[i+1]-rate[i])/rate[i]}

  #Prior distributions for joinpoints' locations
  for(j in 1:4){tau[j]~dunif(tauinf,tausup)}
  tauinf<-x1+2 
  tausup<-x2-2
  #This condition imposes an order restriction and a minimimum distance between joinpoints
  #one has to be given as data and necesarily has to be equal to 1.
  one~dbern(condicion)
  condicion<-step(tau[2]-tau[1]-2)*step(tau[3]-tau[2]-2)*step(tau[4]-tau[3]-2)

  #flat prior distributions for alpha and beta0  
  alpha~dflat()
  beta0~dflat()
  
  #Prior distribution for parameter beta
  beta[1:4]~dmnorm(mean.beta[1:4],prec.u[1:4,1:4])
  for(i in 1:4){
    mean.beta[i]<-0
    for(k in 1:n){BaseJoin2[k,i]<-BaseJoin[k,i]*pow(w[k],1/2)}
    #prior precision matrix for beta.
    #The parameter g makes the prior distribution of beta to be a Cauchy distribution
    for(j in 1:4){
      prec.u[i,j]<-g*inprod(BaseJoin2[,i],BaseJoin2[,j])*max(step(i-j)*step(j-i),delta[i]*delta[j])
    }
  }
  #Prior distribution for g. 
  g~dgamma(0.5,n2)
  n2<-n/2
  
  #Parameter w for the prior precision of beta
  for(j in 1:n){
    w[j]<-Pob[j]*exp(alpha+beta0*(anyo[j]-(x1+x2)/2))
  }
  
  #########################
  #Break-point definitions#
  #########################
  for(j in 1:4){
    for(i in 1:n){
      BaseJoin[i,j]<-1+(anyo[i]-tau[j])*gamma[i,j]
      gamma[i,j]<-gamma1[j]+(gamma2[j]-gamma1[j])*step(anyo[i]-tau[j])
    }
    gamma1[j]<-(d0*c2[j]-c0*d2[j])/(d2[j]*c1[j]-d1[j]*c2[j])
    gamma2[j]<-(d1[j]*c0-d0*c1[j])/(d2[j]*c1[j]-d1[j]*c2[j])
  
    c1[j]<-inprod(anyocent[,j],primerperiodo[,j])
    c2[j]<-inprod(anyocent[,j],segundoperiodo[,j])
    d1[j]<-inprod(anyocent2[,j],primerperiodo[,j])
    d2[j]<-inprod(anyocent2[,j],segundoperiodo[,j])

    for(i in 1:n){
      anyocent[i,j]<-anyo[i]-tau[j]
      anyocent2[i,j]<-anyo[i]*anyocent[i,j]
      primerperiodo[i,j]<-step(-anyocent[i,j])
      segundoperiodo[i,j]<-1-primerperiodo[i,j]
    }

  }
  c0<-n;d0<-sum(anyo[])

  # Prior distributions for delta  
  delta[1]~dbern(p1)  
  delta[2]~dbern(p2)
  delta[3]~dbern(p3)
  delta[4]~dbern(p4) 
  # Prior distributions for the probability of deltas being included in the model
  p1~dbeta(0.5,1.5)
  p2~dbeta(0.5,1.5)
  p3~dbeta(0.5,1.5)
  p4~dbeta(0.5,1.5)

  #Indicator variable pointing which of the models M_0,M_1,M_2,M_3 or M_4 is sampled
  ModelSampled<-sum(delta[])
}