##############################################################################
# Dataset                                                               
##############################################################################

#Dataset 1
D1 = 	c(0,0,0,74.4,0,63.3,36.7,0,0,0
		,0,0,0,10.9,0,0,0,0,4.8,0
		,0,0.3,0,0,0,52.9,7.2,30.5,0,0
		,3.6,20.9,1.9,0,0,0,16.4,0,12,2.8
		,10.4,15.1,0,82,31,29.6,0.9,41.8,0,21.7)
		
#Dataset 2
D2 =	c(0,56.8,66,0,0,0.8,4,1.6,0,0
		,0,1,0,2.2,2.9,0,0,0,0,0
		,0,0,0,13.6,0,1.9,0,47.3,16.5,47.3
		,0,0,24.4,0,6.2,9.2,33.9,21.5,0,60.1
		,32.3,47.3,1,71.5,0,0)

#Dataset 3
D3 =	c(3.8,0,38.3,0,0,0,0,0,4.7,0
		,0,19,12.2,0,21.8,0,23.5,0,0,0
		,31.3,0,0,38.8,32.4,0,4.1,0,0.7,0
		,0,0,8.3,0,0,0,0.5,15.3,0,0
		,0.6,0,0,0,0,50.2)
		
#Dataset 4
D4 =	c(0,0,0,22.5,0,0,3.5,6,13.2,0
		,51.1,29.9,34.3,0,42.5,0,2.5,0,0,0
		,46.8,0,0,0,0,28.7,0,0,7.3,0
		,0,0,0,16,4.7,11,0,5.2,24.5,61.2
		,61.5,7.4,34,0,102.2,0,0,57.2,0,0,0)

#Dataset 5
D5 =	c(0.4,0.1,0.3,1.7,3.3,0,0,0,0,0
		,88.6,0,59.2,34.5,0,0,0,0,0,0
		,7.4,0,0,0,121.5,2.8,0,0,0.5,0
		,0,0,0,0,0,0,7.5,0,25.6,2
		,0,0,0,11,0,0,7.7,12.1,38,0
		,55.5,28,30.4,0,48.4,0,6,13,8.3,0
		,0,0,0,0,0.1,0,0,0,11.4,1.6
		,54.5,0,0,0,0,0,38.3,0,0,0
		,0,0,13.6,0,0,0,0,0,31.7,0
		,0,39,13.2,0,0,0,0.5,0,0,0
		,11,0,0,0,0,27.1,0,0,3.9,0
		,0,0,0,7.4,13.4,0,3.2,1.1,0,43
		,0,12.6,40.4,19.2,0,5.8,29.8,0,0,9.1
		,0,0,0,0.5,35,0,0,0,2.7,9.1
		,0,0,0.8,136.2,0.5,21.6,38.5,0,19.5,44.5
		,83.5,0,24.5,0,29,33.2,1.5,0,26.3,20.2
		,61.4,11.8,71,0,29.1,0,0,15.6,0,0,0.9)
		

##############################################################################
# R code to estimate AIC results                                             
##############################################################################

library(extraDistr)
library(MASS)
library(fitdistrplus)
library(LaplacesDemon)

#Table 4

#Dataset 1
x 	<- D1[which(D1!=0)] 
n1 	<- length(x)
n	<- length(D1)
n0 	<- n-n1
  
delta 	<- n0/(n0+n1)

fit.norm = fitdistr(x,"normal")
fit.lognorm = fitdistr(x,"lognormal")
fit.cauchy = fitdistr(x,"cauchy")
fit.gamma = fitdistr(x,"gamma",lower = c(0, 0),start = list(scale = 1, shape = 1) )  

c(AIC(fit.norm),AIC(fit.lognorm),AIC(fit.cauchy),AIC(fit.gamma))

##############################################################################

#Dataset 2
x 	<- D2[which(D2!=0)] 
n1 	<- length(x)
n	<- length(D2)
n0 	<- n-n1
  
delta 	<- n0/(n0+n1)

fit.norm = fitdistr(x,"normal")
fit.lognorm = fitdistr(x,"lognormal")
fit.cauchy = fitdistr(x,"cauchy")
fit.gamma = fitdistr(x,"gamma",lower = c(0, 0),start = list(scale = 1, shape = 1) )  

c(AIC(fit.norm),AIC(fit.lognorm),AIC(fit.cauchy),AIC(fit.gamma))

##############################################################################

#Dataset 3
x 	<- D3[which(D3!=0)] 
n1 	<- length(x)
n	<- length(D3)
n0 	<- n-n1
  
delta 	<- n0/(n0+n1)

fit.norm = fitdistr(x,"normal")
fit.lognorm = fitdistr(x,"lognormal")
fit.cauchy = fitdistr(x,"cauchy")
fit.gamma = fitdistr(x,"gamma",lower = c(0, 0),start = list(scale = 1, shape = 1) )  

c(AIC(fit.norm),AIC(fit.lognorm),AIC(fit.cauchy),AIC(fit.gamma))

##############################################################################

#Dataset 4
x 	<- D4[which(D4!=0)] 
n1 	<- length(x)
n	<- length(D4)
n0 	<- n-n1
  
delta 	<- n0/(n0+n1)

fit.norm = fitdistr(x,"normal")
fit.lognorm = fitdistr(x,"lognormal")
fit.cauchy = fitdistr(x,"cauchy")
fit.gamma = fitdistr(x,"gamma",lower = c(0, 0),start = list(scale = 1, shape = 1) )  

c(AIC(fit.norm),AIC(fit.lognorm),AIC(fit.cauchy),AIC(fit.gamma))

##############################################################################

#Dataset 5
x 	<- D5[which(D5!=0)] 
n1 	<- length(x)
n	<- length(D5)
n0 	<- n-n1
  
delta 	<- n0/(n0+n1)

fit.norm = fitdistr(x,"normal")
fit.lognorm = fitdistr(x,"lognormal")
fit.cauchy = fitdistr(x,"cauchy")
fit.gamma = fitdistr(x,"gamma",lower = c(0, 0),start = list(scale = 1, shape = 1) )  

c(AIC(fit.norm),AIC(fit.lognorm),AIC(fit.cauchy),AIC(fit.gamma))

##############################################################################