#Define the working directory
setwd("C:/Users/gsyla/OneDrive/Notes/Xanthi")


#Install required packages
install.packages("e1071")
install.packages("fitdistrplus")

#Call package libraries
library(e1071)
library(fitdistrplus)

#read the external data file
flow1 = read.table("Strymon_discharge.txt", 
                   header=TRUE, 
                   dec=".")

#transform list to vector
flow = flow1$Discharge

#Compute descriptive statistics
mean.flow = mean(flow)
median.flow = median(flow)
range = range(flow)
range.flow = range[2]-range[1]
sd.flow = sd(flow)
var.flow = var(flow)
IQR.flow = IQR(flow)
summary(flow)

#Compute skewness and kurtosis
skewness.flow = skewness(flow,type=1)
kurtosis.flow = kurtosis(flow,type=1)

#Plot histogram of Strymon data
hist(flow,
     main="Histogram of Strymon River Discharge", 
     xlab="Discharge, m^3/s", 
     xlim=c(0,600),
     breaks = seq(0,600, by=50),
     col="red")

#Plot boxplot of Strymon data
boxplot(flow,
        main = "Boxplot for Strymon Discharge data",
        xlab = "Strymon",
        ylab = "River Discharge, m^3/s")


#Plot histogram, probability density function
#and cumulative density distribution 
plotdist(flow, histo = TRUE, demp = TRUE)

#Plot the location of Strymon data in relation to other 
#theoretical probability density distributions
descdist(flow)

#Normalize Strymon River data
flow_scaled <- (flow1$Discharge - min(flow1$Discharge) + 0.001) / (max(flow1$Discharge) - min(flow1$Discharge) + 0.002)
summary(flow_scaled)

#Fit the normal probability density distribution 
#to the scaled Strymon River data
fn <- fitdist(flow_scaled,
               "norm",
               method = "mle")

#Get summary of normal model output
summary(fn)

#Make a plot of data and normal distribution curve
plot(fn)

#Fit the gamma probability density distribution
#to the raw Strymon River data
fg <- fitdist(flow_scaled,
              "gamma",
              method = "mle")


summary(fg)
plot(fg)


#Fit the most appropriate Weibull probability density
#distribution model to Strymon data
fw <- fitdist(flow_scaled,
              "weibull",
              method = "mle")


#Get summary of Weibull model output
summary(fw)
plot(fw)

#Denormalize
flow_test = 1000
flow_norm <- (flow_test - min(flow) + 0.001) / (max(flow) - min(flow) + 0.002)


#Get the probability of an event 
#with discharge from 0 to 100 m^3/s with step 50 m^3/s
dgamma(seq(.1,flow_norm,by=0.05), 
         shape=0.573, rate = 4.218,
         log = FALSE)

#Get the cumulative probability of an event 
#with discharge higher than 100 m^3/s 
pgamma(flow_norm, 
       shape=0.573, rate = 4.218, 
       log = FALSE, 
       lower.tail=FALSE)

#Create 10 random values sampled from a population 
#of Weibull probability density distribution
rgamma(n=10, shape=0.573, rate = 4.218)