##### Lesson 55: Continuous Distributions Part II #####

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

# Read the file #
cp_data = read.csv("cp_temperature.csv",header=T)

# Extract March 2017 Daily Temperature Data #
index = which((cp_data[,2]==3) & (cp_data[,1] == 2017))
march_2017 = cp_data[index,4]

library(locfit)
hist(march_2017,prob=T,main="",font=2,font.lab=2,xlab="March 2017 Daily Temperature")
lines(locfit(~march_2017))

# Extract March 2016 Daily Temperature Data #
index = which((cp_data[,2]==3) & (cp_data[,1] == 2016))
march_2016 = cp_data[index,4]

hist(march_2016,prob=T,main="",font=2,font.lab=2,xlab="March 2016 Daily Temperature")
lines(locfit(~march_2016))


## Code to demonstrate the Central Limit Theorem ##
yr1 = 1951
yr2 = 2017
n = yr2 - yr1 + 1

index = which((cp_data[,2]==3) & (cp_data[,1] == yr1)) 
march_daily_temperature = cp_data[index,4] 

par(mfrow=c(2,1))
plot(locfit(~march_daily_temperature),main="Average Daily Temperature",col="grey",xlim=c(10,80),ylim=c(0,0.1),font=2,font.lab=2,xlab="Temperature (deg F)")
points(mean(march_daily_temperature,na.rm=T),0,col="red",pch=25)

avg = matrix(NA,nrow=n,ncol=1)
for (i in 1: n)
{
  yr = 1950 + i
  index2 = which((cp_data[,1]==yr) & (cp_data[,2]==3))
  daily =  cp_data[index2,4]

  lines(locfit(~daily),col="grey")
  
  avg[i,1] = mean(cp_data[index2,4],na.rm=T)
  points(avg[i,1],0,col="red",pch=25)
  Sys.sleep(0.5) 
}
plot(locfit(~avg),main="Average Monthly Temperature",col="white",lwd=2,xlim=c(10,80),font=2,font.lab=2,xlab="Temperature (deg F)")
lines(locfit(~avg),col="red",lwd=2)

## Normal Distribution ##
hist(avg,prob=T,main="Average Monthly Temperature",lwd=2,font=2,font.lab=2,xlab="Temperature (deg F)")
lines(locfit(~avg),col="grey")

# Probability Computations #
pnorm(40,mean(avg),sd(avg))

1 - pnorm(50,mean(avg),sd(avg))

# Quantiles # 
qnorm(0.1,mean(avg),sd(avg))

## Lognormal Distribution #

## Vapor Pressure ##
Tc = (avg-32)*(5/9) 
vp = 0.61078*exp((17.27*Tc)/(Tc+243.05))

# Plotting and Comparing #
par(mfrow=c(2,1))
hist(Tc,prob=T,main="Average Monthly Temperature",lwd=2,font=2,font.lab=2,xlab="Temperature (deg C)")
lines(locfit(~Tc),col="grey")

hist(vp,prob=T,main="Vapor Pressure ~ exp(T)",lwd=2,font=2,font.lab=2,xlab="Vapor Pressure (kpa)")
lines(locfit(~vp),col="grey")