#
# Change-point analysis of coal disaster data
#
m <- 2                           #  Number of replications
t <- 15                          #  Number of iterations
startyear <- 1851                #  First year for which data is available
daytab  <- c(0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
leaptab <- c(0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)

#  Functions day.of.year, month.of.day and this.month adapted from 
#  B W Kernighan and D M Ritchie, The C Programming Language,
#  Englewood Cloffs, NJ: Prentice-Hall 1978, 1988, Section 5.7.

#  day.of.year: set day of year from month & day
day.of.year <- function(year,month,day)
{
  leap <- year%%4 == 0 && year%%100 != 0 || year%%400 == 0
  if (leap)
    tab <- leaptab else
    tab <- daytab
  yearday <- day
  for (i in 1:month)
    yearday <- yearday + tab[i]
  return(yearday)
}

#  month.of.day: set month, day from day of year
day.of.month <- function(year,yearday)
{
  leap <- year%%4 == 0 && year%%100 != 0 || year%%400 == 0
  if (leap)
    tab <- leaptab else
    tab <- daytab
  i <- 1
  while (yearday > tab[i])
    {
      yearday <- yearday - tab[i]
      i <- i + 1
    }
  return(yearday)
}

this.month <- function(year,yearday)
{
  leap <- year%%4 == 0 && year%%100 != 0 || year%%400 == 0
  if (leap)
    tab <- leaptab else
    tab <- daytab
  i <- 1
  while (yearday > tab[i])
    {
      yearday <- yearday - tab[i]
      i <- i + 1
    }
  return(i-1)
}

#  Data from B P Carlin, A E Gelfand and A F M Smith, Hierachical 
#  Bayesian Analysis of Changepoint Problems, Appl. Statist. 41 (1992), 
#  389-405.
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)                   #  Number of years of data available
endyear <- startyear+n-1         #  First year for which data is available
a1 <- 0.5
a2 <- 0.5
plot(startyear:endyear,cumsum(Y))
cat("\n")
pp <- rep(0,n)
L <- rep(0,n)
pp <- rep(0,n)
for (j in 1:m)                   #  Replicate m times
  {
    k <- 1+floor(n*runif(1))     #  Initialize k randomly in [1,n]
    gamma <- 1
    delta <- 1                      #  Initialize b1=b2=1
    for (s in 1:t)               #  Iterate t times
      {
        #  Sample lambda | Y,mu,gamma,delta,k
        lambda <- rgamma(1,a1+cumsum(Y)[k])/(gamma+k)
        #  Sample mu | Y,lambda,gamma,delta,k
        mu <- rgamma(1,a2+sum(Y)-cumsum(Y)[k])/(delta+n-k)
        #  Sample gamma | Y,lambda,mu,delta,k
        gamma <- rgamma(1,a1)/(1+lambda)
        #  Sample delta | Y,lambda,mu,gamma,k
        delta <- rgamma(1,a2)/(1+mu)
        #  Find L(k | lambda,mu,Y) for k = 0 to n-1
        for (k in 1:n)
          {
             L[k] <- exp(k*(mu-lambda)+
                         (log(lambda)-log(mu))*cumsum(Y)[k])
          }
        #  Find p(k | lambda,mu,gamma,delta) and cumulation thereof
        p <- L/sum(L)
        cumprob <- cumsum(p)
        #  Pick U at random between 0 and 1
        U <- runif(1)
        #  Sample k | Y,lambda,mu,gamma,delta
        for (i in 1:n)
          if ((cumprob[i] < U)&&(U <= cumprob[i+1])) k <- i
      } # End iteration
  pp <- pp + p/m
  }     # End replication
#  Find posterior density and mean of k
year <- startyear:endyear
meandate <- sum((year+0.5)*pp)
#  Print out results
for (i in 30:50) cat(startyear+i," ",pp[i],"\n")
cat("\n")
for (i in 30:50)
  {
    cat(startyear+i," ")
    for (j in 1:80)
      if (100*pp[i] > j) cat("*")
    cat("\n")
  }
cat("\n")
meanyear <- floor(meandate)
fracyear <- meandate - floor(meandate)
leap <- meanyear%%4 == 0 && meanyear%%100 != 0 || meanyear%%400 == 0
if (leap)
  tab <- leaptab else
  tab <- daytab
daysinyear <- if (leap) 366 else 365
remnant <- fracyear*(daysinyear)-cumsum(tab)
monthspast <- remnant[remnant>0]
meanmonth <- length(monthspast)
if (meanmonth==1) monthname <- "Jan"
if (meanmonth==2) monthname <- "Feb"
if (meanmonth==3) monthname <- "Mar"
if (meanmonth==4) monthname <- "Apr"
if (meanmonth==5) monthname <- "May"
if (meanmonth==6) monthname <- "Jun"
if (meanmonth==7) monthname <- "Jul"
if (meanmonth==8) monthname <- "Aug"
if (meanmonth==9) monthname <- "Sep"
if (meanmonth==10) monthname <- "Oct"
if (meanmonth==11) monthname <- "Nov"
if (meanmonth==12) monthname <- "Dec"
floatday <- fracyear*(daysinyear)-cumsum(tab)[meanmonth]
meanday <- 1+floor(floatday)
cat("Mean is",meanday,monthname,meanyear,", i.e. ")
cat(meanyear,"+",fracyear,"\n")
cat("\n")