model{
  mu ~ dnorm(0,0.001)
  tau ~dgamma(0.001,0.001)
  for(i in 1:n){
    y[i] ~ dnorm(mu,tau)
    yrep[i] ~ dnorm(mu,tau)
  }
  n.50 <- round(n/2)
  n.25 <- round(n/4)
  yrep.sort <- sort(yrep[])
  T1.rep <- yrep.sort[1]
  yrep.50 <- yrep.sort[n.50]
  yrep.25 <- yrep.sort[n.25]
  T2.rep <- (T1.rep-yrep.50)/(yrep.25-yrep.50)
  P1 <- step(T1.rep-T1.obs)
  P2 <- step(T2.rep-T2.obs)
}