R Under development (unstable) (2024-07-10 r86888 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(fitdistrplus)
Loading required package: MASS
Loading required package: survival
> 
> set.seed(123)
> 
> nsample <- 10
> 
> #--------------------------------------------------------
> # exponential sample
> 
> x1 <- rexp(nsample)
> 
> mseKL_exp <- msedist(x1, "exp")
> mseJ_exp <- msedist(x1, "exp", phidiv = "J")
> mseR2_exp <- msedist(x1, "exp", phidiv = "R", power.phidiv=2)
> mseR1o2_exp <- msedist(x1, "exp", phidiv = "R", power.phidiv=1/2)
> mseH3o2_exp <- msedist(x1, "exp", phidiv = "H", power.phidiv=3/2)
> mseV3o2_exp <- msedist(x1, "exp", phidiv = "V", power.phidiv=3/2)
> 
> c(true=1, mseKL_exp$estimate, mseJ_exp$estimate, mseR2_exp$estimate, 
+   mseR1o2_exp$estimate, mseH3o2_exp$estimate, mseV3o2_exp$estimate)
    true     rate     rate     rate     rate     rate     rate 
1.000000 1.336924 1.329121 1.164714 1.195323 1.167898 1.171918 
> 
> 
> 
> mseKL_exp <- fitdist(x1, "exp", method="mse", phidiv="KL")
> mseJ_exp <- fitdist(x1, "exp", method="mse", phidiv = "J")
> mseR2_exp <- fitdist(x1, "exp", method="mse", phidiv = "R", power.phidiv=2)
> mseR1o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "R", power.phidiv=1/2)
> mseH3o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "H", power.phidiv=3/2)
> mseV3o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "V", power.phidiv=3/2)
> 
> 
> 
> gofstat(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp))
Goodness-of-fit statistics
                             1-mse-exp 2-mse-exp 3-mse-exp 4-mse-exp 5-mse-exp
Kolmogorov-Smirnov statistic 0.2549881 0.2566077 0.2916788 0.2850105 0.2909822
Cramer-von Mises statistic   0.1544480 0.1568271 0.2185483 0.2052244 0.2171190
Anderson-Darling statistic   1.0571407 1.0703179 1.4200650 1.3435948 1.4118434
                             6-mse-exp
Kolmogorov-Smirnov statistic 0.2901037
Cramer-von Mises statistic   0.2153293
Anderson-Darling statistic   1.4015539

Goodness-of-fit criteria
                               1-mse-exp 2-mse-exp 3-mse-exp 4-mse-exp
Akaike's Information Criterion  13.22043  13.23811  13.78497  13.65601
Bayesian Information Criterion  13.52301  13.54070  14.08756  13.95860
                               5-mse-exp 6-mse-exp
Akaike's Information Criterion  13.77093  13.75341
Bayesian Information Criterion  14.07351  14.05599
> 
> 
> cdfcomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), do.points=FALSE,
+         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
+                        "Vajda beta=3/2"))
> 
> qqcomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), 
+         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
+                        "Vajda beta=3/2"))
> 
> denscomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), demp = TRUE,
+        legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
+                       "Vajda beta=3/2"))
> 
> 
> #defensive test 
> 
> try(msedist(x1, "exp", phidiv="ABC"))
Error in match.arg(phidiv, c("KL", "J", "R", "H", "V")) : 
  'arg' should be one of "KL", "J", "R", "H", "V"
> try(msedist(x1, "exp", phidiv="K", power.phidiv = "a"))
Error in msedist(x1, "exp", phidiv = "K", power.phidiv = "a") : 
  is.null(power.phidiv) is not TRUE
> try(msedist(x1, "exp", phidiv="J", power.phidiv = "a"))
Error in msedist(x1, "exp", phidiv = "J", power.phidiv = "a") : 
  is.null(power.phidiv) is not TRUE
> try(msedist(x1, "exp", phidiv="R", power.phidiv = 0))
Error in msedist(x1, "exp", phidiv = "R", power.phidiv = 0) : 
  power.phidiv > 0 && power.phidiv != 1 is not TRUE
> try(msedist(x1, "exp", phidiv="R", power.phidiv = 1:10))
Error in msedist(x1, "exp", phidiv = "R", power.phidiv = 1:10) : 
  length(power.phidiv) == 1 is not TRUE
> try(msedist(x1, "exp", phidiv="H", power.phidiv = 0))
Error in msedist(x1, "exp", phidiv = "H", power.phidiv = 0) : 
  power.phidiv >= 1 is not TRUE
> try(msedist(x1, "exp", phidiv="H", power.phidiv = 1:10))
Error in msedist(x1, "exp", phidiv = "H", power.phidiv = 1:10) : 
  length(power.phidiv) == 1 is not TRUE
> 
> 
> 
> #--------------------------------------------------------
> # Poisson sample
> 
> x1 <- rpois(nsample, lambda=15)
> 
> #no weight
> mseKL_pois1 <- fitdist(x1, "pois", method="mse", phidiv="KL")
> mseJ_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "J")
> mseR2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "R", power.phidiv=2)
> mseR1o2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "R", power.phidiv=1/2)
> mseH3o2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "H", power.phidiv=3/2)
> mseV3o2_pois1 <- fitdist(x1, "pois",method="mse",  phidiv = "V", power.phidiv=3/2)
> 
> 
> #with weight
> mseKL_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv="KL", weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> mseJ_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "J", weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> mseR2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "R", power.phidiv=2, weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> mseR1o2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "R", power.phidiv=1/2, weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> mseH3o2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "H", power.phidiv=3/2, weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> mseV3o2_pois2 <- fitdist(unique(sort(x1)), "pois",method="mse",  phidiv = "V", power.phidiv=3/2, weights=as.numeric(table(x1)))
Warning message:
In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
  weights are not taken into account in the default initial values
> 
> #taking into account partially unbias the estimation of the true cdf
> cdfcomp(list(mseKL_pois1, mseJ_pois1, mseR2_pois1), do.points=FALSE, fitlty=1,
+         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2"))
> 
> cdfcomp(list(mseKL_pois2, mseJ_pois2, mseR2_pois2), do.points=FALSE, add=TRUE, 
+         fitlty=2, addlegend = FALSE, datacol="white")
> 
> 
> cdfcomp(list(mseR1o2_pois1, mseH3o2_pois1, mseV3o2_pois1), do.points=FALSE, fitlty=1,
+         legendtext = c("Renyi alpha=1/2", "Hellinger p=3/2", "Vajda beta=3/2"))
> 
> cdfcomp(list(mseR1o2_pois2, mseH3o2_pois2, mseV3o2_pois2), do.points=FALSE, add=TRUE, 
+         fitlty=2, addlegend = FALSE, datacol="white")
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.73    0.32    1.98