##### Lesson 6: It’s “probably” me #####

# Setting the working directory 
setwd(“provide the path to your folder”)
list.files()

# reading the input files
jfktemperature = read.csv("jfk_dailytemperature.csv",header=TRUE)

# number of years in the data 
dum = range(jfktemperature[,1]) # extracting the range of years in the data 
years = (dum[1]+1):dum[2] # creating a sequence of years from 1948 to 2017
years = rev(years)

# for each year count the number of warm days (temperature > 50 F or 10 C) in February

# creating a matrix with rows = number of years and 5 columns; first column = year 
# second column = days in february, third column = warm days in that year, 
# fourth column = cumulative days in february, fifth column = cumulative warm days 
# sixth column = probability of warm days 
warm_threshold = 10
nyrs = length(years)
feb_warmdays = matrix(NA,nrow=nyrs,ncol=6)
for (i in 1:nyrs)
{
  dum = which( (jfktemperature[,1]==years[i]) & (jfktemperature[,2]==2))
  warmdays = length(which(jfktemperature[dum,4]>=warm_threshold))
  feb_warmdays[i,1] = years[i]
  feb_warmdays[i,2] = length(dum)
  feb_warmdays[i,3] = (warmdays)
  feb_warmdays[i,4] = sum(feb_warmdays[1:i,2])
  feb_warmdays[i,5] = sum(feb_warmdays[1:i,3])
  feb_warmdays[i,6] = feb_warmdays[i,5]/feb_warmdays[i,4]
}
  
plot(feb_warmdays[,6],type="o",lty=2,lwd=2,axes=F,font.lab=2,xlab="Years upto",ylab="Probability of Warm Days in February")
axis(1, at=1:nyrs, years,las = 1, cex.axis = 1,font=2); 
axis(2,at=seq(0,0.25,by=0.02),las=1, cex.axis = 1,font=2)
abline(h=0.027,lwd=2,lty=2,col="red")