# Supplementary material to the manuscript
# "Calculating nuclear accidents probabilities from empirical frequencies"
# Minh Ha-Duong, CNRS      and  Venance Journé, CNRS
# haduong@centre-cired.fr       journe@centre-cired.fr
#
# (c) Ha-Duong, 2012, 2013
# This R script and the accompaning data file are released according to Creative Commons "Attribution 3.0" license

yearStart <- 1951
yearINES  <- 1991
yearNow   <- 2011
year2date <- function(year) {(year-1970)*365.25}

data <- read.csv("accidents.csv" , header = TRUE, sep="\t")
attach(data)

#Table 2: Number of events by INES level

table(data$INES, useNA="always")

table(subset(data,FACILITY == "Power reactor")$INES, useNA="always")
table(subset(data,FACILITY != "Power reactor")$INES, useNA="always")

table(subset(data, as.Date(DATE) >= year2date(yearINES))$INES, useNA="always")
table(subset(data, as.Date(DATE) < year2date(yearINES))$INES, useNA="always")


# Table 3  Frequency distributions of the number of occurences during a set of random years

# Run once if needed, as root if needed
#install.packages("polynom")

require(polynom)

probaTable <- function(observations, sampleSize, minimalGravity, maxNi) {
  events.count <- tabulate(subset(observations$YEAR, observations$INES>=minimalGravity))[yearStart:yearNow] ;
  events.table <- table(factor(events.count, levels=c(0,1,2,3))) / (yearNow-yearStart+1) ;
  events.generatingFunction  <- polynomial(events.table) ;
  sampleSum.generatingFunction <- events.generatingFunction ** sampleSize;
  sampleSum.probas = as.vector(sampleSum.generatingFunction);
  result   <- c(sampleSum.probas[1:maxNi], sum(sampleSum.probas[(maxNi+1):length(sampleSum.probas)])) ;
  return(round(result,3)*100)
}

reactorData=subset(data,FACILITY=="Power reactor")
fukushimaEvent=subset(reactorData, LOCATION=="Fukushima daichi")
sensitivityData=rbind(rbind(reactorData, fukushimaEvent),fukushimaEvent)

probaTable(reactorData, 5, 4, 4)
probaTable(reactorData, 5, 7, 4)
probaTable(reactorData, 30, 4, 9)
probaTable(reactorData, 30, 7, 4)
probaTable(sensitivityData, 5, 4, 4)
probaTable(sensitivityData, 5, 7, 4)


#Double checking Table 3 with a different method

require(boot); nReplicates <- 100000

probaTable <- function(observations, sampleSize, minimalGravity, maxNi) {
  headSum <- function(d,i) sum(d[i[1:sampleSize]])
  events.count <- tabulate(subset(observations$YEAR, observations$INES>=minimalGravity))[yearStart:yearNow] ;
  events.boot  <- boot(data=events.count, statistic=headSum, R=nReplicates) ;
  events.table <- table(as.vector(events.boot[2])) ;
  events.dist  <- events.table / nReplicates ;
  events.row   <- c(events.dist[1:maxNi], sum(events.dist[(maxNi+1):length(events.dist)])) ;
  names(events.row)[maxNi+1] <- paste(maxNi,'+');
  return(round(events.row,3)*100)
}

probaTable(reactorData, 5, 4, 4)
probaTable(reactorData, 5, 7, 4)
probaTable(reactorData, 30, 4, 9)
probaTable(reactorData, 30, 7, 4)
probaTable(sensitivityData, 5, 4, 4)
probaTable(sensitivityData, 5, 7, 4)