set.seed(1234)
library(MASS)
library(vcd)
library(survival)
library(fitdistrplus)
library(logspline)
library(actuar)
library(fitur)
library(distr)
library(tolerance)
library(Renext)
library(nsRFA)

#===Data Rainfall Kiew Lom Dam January 1993-2021===#
x1 = c(0,0.1,0,0,0,0,0,0,15.7,3.7,0,0,0,0,0,34.5,0,43.4,4.3,34.6,38.5,0,29.5,65.7,51.9,21.7,14.9,0,0)

#===Data Rainfall Kiew Lom Dam February 1993-2021===#
x2 = c(0,0,22.0,2,0,0,0,20.2,0,0.7,0,0,13,9.2,0,11.2,0,0,0,0.5,4.5,0,0,3.5,0,19.5,12,0,8)

#===Data Rainfall Kiew Lom Dam January-February 1993-2021===#
x = c(x1,x2)

#----------------Shape, Rate, scale=1/rate----------------------------------------------------------------#
n = length(x)
nonzero.value = x[which(x!=0)]
x.nonzero <- x[which(x!=0)]
n.1 = length(nonzero.value)
n.0 = n-n.1
n
n.1
n.0
x.bar = mean(nonzero.value)
x.bar
shape.hat =  1/(2*(log(x.bar)-(sum(log(nonzero.value))/n.1)))
rate.hat = shape.hat/x.bar
scale.hat = 1/rate.hat
delta.hat = n.0/n
shape.hat
rate.hat
scale.hat
delta.hat

meanDG.hat <- (1-delta.hat)*(shape.hat/rate.hat)
tau.hat <- ((1-delta.hat)*(shape.hat/(rate.hat^2)))+((delta.hat*(1-delta.hat))*((shape.hat/rate.hat)^2))
alpha <- 0.05
meanDG.hat
tau.hat
#--------------------------------------------------------------------------------------------------------#
str(x)
plotdist(x)
plotdist(x, histo = TRUE, demp = TRUE)
descdist(x, boot=1000)
descdist(x,discrete=FALSE)
qqcomp(fitgamma)

#------Fit distribution------#
fitnormal <- fitdist(nonzero.value, "norm")
fitcauchy <- fitdist(nonzero.value, "cauchy")
fitlognormal <- fitdist(nonzero.value, "lnorm")
fitgamma <- fitdist(nonzero.value, "gamma", lower=c(0,0), start=list(scale=1, shape=1))

summary(fitnormal)
summary(fitcauchy)
summary(fitlognormal)
summary(fitgamma)

plot(fitgamma)