##### Lesson 55: Continuous Distributions Part I #####

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

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

# Select Wind Gust data > 35 mph from the file #
gust_events = length(which(cp_data$WSFG > 35))

# Count the number of events each year #
years = range(cp_data$Year)

number_years = diff(years) + 1

gust_counts_yearly = matrix(NA,nrow=number_years,ncol=1)

for (i in 1:number_years)
{
  year = 1983 + i
  index = which(cp_data$Year == year)
  gust_counts_yearly[i,1] = length(which(cp_data$WSFG[index] > 35))
}

plot(1984:2017,gust_counts_yearly,type="o",font=2,font.lab=2,xlab="Years",ylab="Gust Counts")

hist(gust_counts_yearly)
lambda = mean(gust_counts_yearly)

# create one scenario #
counts_scenario = rpois(34,lambda)


### Beta / Uniform Distribution ###

### Animation ###

a = seq(from = 0.1, to = 2, by = 0.1)
b = 1

x = seq(from = 0, to = 1, by = 0.01)
for (i in 1:length(a))
{
  c = factorial(a[i]+b-1)/(factorial(a[i]-1)*factorial(b-1))
  fx = c*x^(a[i]-1)*(1-x)^(b-1)
  
  plot(x,fx,type="l",font=2,font.lab=2,ylab="f(x)",ylim=c(0,5),main="Beta Family")
  txt1 = paste("a = ",a[i])
  txt2 = paste("b = ",b)
  
  text(0.6,4,txt1,col="red",font=2)
  text(0.6,3.5,txt2,col="red",font=2)
  
  Sys.sleep(1) 
}

a = 2
b = seq(from = 1, to = 2, by = 0.1)

x = seq(from = 0, to = 1, by = 0.01)
for (i in 1:length(b))
{
  c = factorial(a+b[i]-1)/(factorial(a-1)*factorial(b[i]-1))
  fx = c*x^(a-1)*(1-x)^(b[i]-1)
  
  plot(x,fx,type="l",font=2,font.lab=2,ylab="f(x)",ylim=c(0,5),main="Beta Family")
  txt1 = paste("a = ",a)
  txt2 = paste("b = ",b[i])
  
  text(0.6,4,txt1,col="red",font=2)
  text(0.6,3.5,txt2,col="red",font=2)
  
  Sys.sleep(1) 
}

### Exponential Distribution ###

# Select the indices of the Wind Gust data > 35 mph from the file #
library(locfit)

gust_index = which(cp_data$WSFG > 35)

gust_waittime = diff(gust_index)

x = as.numeric(gust_waittime)
hist(x,prob=T,main="",xlim=c(0.1,200),lwd=2,xlab="Wait time in days",ylab="f(x)",font=2,font.lab=2)
lines(locfit(~x),xlim=c(0.1,200),col="red",lwd=2,xlab="Wait time in days",ylab="f(x)",font=2,font.lab=2)

ave_waittime = mean(gust_waittime)

# Simulating wait times #
simulated_waittimes = rexp(34,lambda)*365
plot(locfit(~simulated_waittimes),xlim=c(0.1,200),col="red",lwd=2,xlab="Wait time in days",ylab="f(x)",font=2,font.lab=2)

# Probability of wait time greater than 20 days #
lambda = 1/ave_waittime

1 - pexp(20, lambda)

## Memoryless property ##

# animation plot to understand the memoryless property #
Sys.sleep(1)
t = seq(from=0,to=200,by=1)
ft = lambda*exp(-lambda*t)
plot(t,ft,type="l",font=2,font.lab=2,xlab="Wait Time to Gusty Events (days)",ylab="f(t)")

t = seq(from=20,to=200,by=1)

Sys.sleep(1) 

abline(v=ave_waittime,col="red")

Sys.sleep(1) 

ft = lambda*exp(-lambda*t)
lines(t,ft,col="red",lty=2,lwd=2,type="o")

Sys.sleep(1) 

t = seq(from = 20, to = 200, by =1)
ft = lambda*exp(-lambda*(t-20))
lines(t,ft,col="red",type="o")

#### Gamma Distribution ###

r = seq(from=1,to=5,by=1)
t = seq(from=0,to=200,by=1)

for (i in 1:length(r))
{
  gt = (1/factorial(r[i]-1))*(lambda*exp(-lambda*t)*(lambda*t)^(r[i]-1))
  
  plot(t,gt,type="l",font=2,font.lab=2,xlab="Wait time",ylab="f(t)")
  txt1 = paste("lambda = ",round(lambda,2),sep="")
  txt2 = paste("r = ",r[i],sep="")
  text(50,0.002,txt1,col="red",font=2,font.lab=2)
  text(50,0.004,txt2,col="red",font=2,font.lab=2)
  Sys.sleep(1) 
}

# alternate way #
x1 = rgamma(10000,shape=1,scale=ave_waittime)
plot(locfit(~x1),font=2,font.lab=2,xlab="Wait time",ylab="f(t)")
txt1 = paste("lambda = ",round(lambda,2),sep="")
txt2 = paste("r = ",1,sep="")
text(50,0.002,txt1,col="red",font=2,font.lab=2)
text(50,0.004,txt2,col="red",font=2,font.lab=2)

x2 = rgamma(10000,shape=2,scale=ave_waittime)
plot(locfit(~x2),font=2,font.lab=2,xlab="Wait time",ylab="f(t)")
txt1 = paste("lambda = ",round(lambda,2),sep="")
txt2 = paste("r = ",2,sep="")
text(50,0.002,txt1,col="red",font=2,font.lab=2)
text(50,0.004,txt2,col="red",font=2,font.lab=2)

x3 = rgamma(10000,shape=3,scale=ave_waittime)
plot(locfit(~x3),font=2,font.lab=2,xlab="Wait time",ylab="f(t)")
txt1 = paste("lambda = ",round(lambda,2),sep="")
txt2 = paste("r = ",3,sep="")
text(50,0.002,txt1,col="red",font=2,font.lab=2)
text(50,0.004,txt2,col="red",font=2,font.lab=2)