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
> 
> # ?denscomp
> visualize <- FALSE # TRUE for manual test with visualization of plots
> nsample <- 1000
> nsample <- 10
> 
> # (1) Plot various distributions fitted to serving size data
> #
> data(groundbeef)
> serving <- groundbeef$serving
> fitW <- fitdist(serving,"weibull")
> fitln <- fitdist(serving,"lnorm")
> fitg <- fitdist(serving,"gamma")
> 
> #sanity checks
> try(denscomp("list(fitW, fitln, fitg)",horizontals = FALSE), silent=TRUE)
> try(denscomp(list(fitW, fitln, fitg, a=1),horizontals = FALSE), silent=TRUE)
> 
> #real call
> res <- denscomp(list(fitW, fitln, fitg), probability = TRUE)
> str(res)
List of 2
 $ hist     :List of 6
  ..$ breaks  : num [1:11] 0 20 40 60 80 100 120 140 160 180 ...
  ..$ counts  : int [1:10] 11 37 68 67 20 17 27 4 0 3
  ..$ density : num [1:10] 0.00217 0.00728 0.01339 0.01319 0.00394 ...
  ..$ mids    : num [1:10] 10 30 50 70 90 110 130 150 170 190
  ..$ xname   : chr "mydata"
  ..$ equidist: logi TRUE
  ..- attr(*, "class")= chr "histogram"
 $ densities: num [1:101, 1:3] 0.0021 0.00257 0.00305 0.00353 0.00401 ...
> 
> denscomp(list(fitW, fitln, fitg), probability = FALSE)
> if (requireNamespace ("ggplot2", quietly = TRUE)) {
+   
+   res <- denscomp(list(fitW, fitln, fitg), probability=TRUE, plotstyle = "ggplot")
+   str(res)
+   
+   
+   denscomp(list(fitW, fitln, fitg), probability=FALSE, plotstyle = "ggplot")
+ }
List of 11
 $ data       :'data.frame':	303 obs. of  3 variables:
  ..$ values: num [1:303] 0.0021 0.00257 0.00305 0.00353 0.00401 ...
  ..$ ind   : Factor w/ 3 levels "weibull","lnorm",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..$ sfin  : num [1:303] 10 11.9 13.8 15.7 17.6 19.5 21.4 23.3 25.2 27.1 ...
 $ layers     :List of 2
  ..$ :Classes 'LayerInstance', 'Layer', 'ggproto', 'gg' <ggproto object: Class LayerInstance, Layer, gg>
    aes_params: list
    compute_aesthetics: function
    compute_geom_1: function
    compute_geom_2: function
    compute_position: function
    compute_statistic: function
    computed_geom_params: NULL
    computed_mapping: NULL
    computed_stat_params: NULL
    constructor: call
    data: data.frame
    draw_geom: function
    finish_statistics: function
    geom: <ggproto object: Class GeomBar, GeomRect, Geom, gg>
        aesthetics: function
        default_aes: uneval
        draw_group: function
        draw_key: function
        draw_layer: function
        draw_panel: function
        extra_params: just na.rm orientation
        handle_na: function
        non_missing_aes: xmin xmax ymin ymax
        optional_aes: 
        parameters: function
        rename_size: TRUE
        required_aes: x y
        setup_data: function
        setup_params: function
        use_defaults: function
        super:  <ggproto object: Class GeomRect, Geom, gg>
    geom_params: list
    inherit.aes: TRUE
    layer_data: function
    map_statistic: function
    mapping: uneval
    position: <ggproto object: Class PositionStack, Position, gg>
        compute_layer: function
        compute_panel: function
        fill: FALSE
        required_aes: 
        reverse: FALSE
        setup_data: function
        setup_params: function
        type: NULL
        vjust: 1
        super:  <ggproto object: Class Position, gg>
    print: function
    setup_layer: function
    show.legend: FALSE
    stat: <ggproto object: Class StatBin, Stat, gg>
        aesthetics: function
        compute_group: function
        compute_layer: function
        compute_panel: function
        default_aes: uneval
        dropped_aes: weight
        extra_params: na.rm orientation
        finish_layer: function
        non_missing_aes: 
        optional_aes: 
        parameters: function
        required_aes: x|y
        retransform: TRUE
        setup_data: function
        setup_params: function
        super:  <ggproto object: Class Stat, gg>
    stat_params: list
    super:  <ggproto object: Class Layer, gg> 
  ..$ :Classes 'LayerInstance', 'Layer', 'ggproto', 'gg' <ggproto object: Class LayerInstance, Layer, gg>
    aes_params: list
    compute_aesthetics: function
    compute_geom_1: function
    compute_geom_2: function
    compute_position: function
    compute_statistic: function
    computed_geom_params: NULL
    computed_mapping: NULL
    computed_stat_params: NULL
    constructor: call
    data: data.frame
    draw_geom: function
    finish_statistics: function
    geom: <ggproto object: Class GeomLine, GeomPath, Geom, gg>
        aesthetics: function
        default_aes: uneval
        draw_group: function
        draw_key: function
        draw_layer: function
        draw_panel: function
        extra_params: na.rm orientation
        handle_na: function
        non_missing_aes: linewidth colour linetype
        optional_aes: 
        parameters: function
        rename_size: TRUE
        required_aes: x y
        setup_data: function
        setup_params: function
        use_defaults: function
        super:  <ggproto object: Class GeomPath, Geom, gg>
    geom_params: list
    inherit.aes: TRUE
    layer_data: function
    map_statistic: function
    mapping: uneval
    position: <ggproto object: Class PositionIdentity, Position, gg>
        compute_layer: function
        compute_panel: function
        required_aes: 
        setup_data: function
        setup_params: function
        super:  <ggproto object: Class Position, gg>
    print: function
    setup_layer: function
    show.legend: NA
    stat: <ggproto object: Class StatIdentity, Stat, gg>
        aesthetics: function
        compute_group: function
        compute_layer: function
        compute_panel: function
        default_aes: uneval
        dropped_aes: 
        extra_params: na.rm
        finish_layer: function
        non_missing_aes: 
        optional_aes: 
        parameters: function
        required_aes: 
        retransform: TRUE
        setup_data: function
        setup_params: function
        super:  <ggproto object: Class Stat, gg>
    stat_params: list
    super:  <ggproto object: Class Layer, gg> 
 $ scales     :Classes 'ScalesList', 'ggproto', 'gg' <ggproto object: Class ScalesList, gg>
    add: function
    add_defaults: function
    add_missing: function
    backtransform_df: function
    clone: function
    find: function
    get_scales: function
    has_scale: function
    input: function
    map_df: function
    n: function
    non_position_scales: function
    scales: list
    train_df: function
    transform_df: function
    super:  <ggproto object: Class ScalesList, gg> 
 $ guides     :Classes 'Guides', 'ggproto', 'gg' <ggproto object: Class Guides, gg>
    add: function
    assemble: function
    build: function
    draw: function
    get_custom: function
    get_guide: function
    get_params: function
    get_position: function
    guides: list
    merge: function
    missing: <ggproto object: Class GuideNone, Guide, gg>
        add_title: function
        arrange_layout: function
        assemble_drawing: function
        available_aes: any
        build_decor: function
        build_labels: function
        build_ticks: function
        build_title: function
        draw: function
        draw_early_exit: function
        elements: list
        extract_decor: function
        extract_key: function
        extract_params: function
        get_layer_key: function
        hashables: list
        measure_grobs: function
        merge: function
        override_elements: function
        params: list
        process_layers: function
        setup_elements: function
        setup_params: function
        train: function
        transform: function
        super:  <ggproto object: Class GuideNone, Guide, gg>
    package_box: function
    print: function
    process_layers: function
    setup: function
    subset_guides: function
    train: function
    update_params: function
    super:  <ggproto object: Class Guides, gg> 
 $ mapping    :List of 4
  ..$ x     : language ~.data$sfin
  .. ..- attr(*, ".Environment")=<environment: 0x0000025b30e02ea0> 
  ..$ y     : language ~.data$values
  .. ..- attr(*, ".Environment")=<environment: 0x0000025b30e02ea0> 
  ..$ group : language ~.data$ind
  .. ..- attr(*, ".Environment")=<environment: 0x0000025b30e02ea0> 
  ..$ colour: language ~.data$ind
  .. ..- attr(*, ".Environment")=<environment: 0x0000025b30e02ea0> 
  ..- attr(*, "class")= chr "uneval"
 $ theme      :List of 2
  ..$ plot.title     :List of 11
  .. ..$ family       : NULL
  .. ..$ face         : NULL
  .. ..$ colour       : NULL
  .. ..$ size         : NULL
  .. ..$ hjust        : num 0.5
  .. ..$ vjust        : NULL
  .. ..$ angle        : NULL
  .. ..$ lineheight   : NULL
  .. ..$ margin       : NULL
  .. ..$ debug        : NULL
  .. ..$ inherit.blank: logi FALSE
  .. ..- attr(*, "class")= chr [1:2] "element_text" "element"
  ..$ legend.position: chr "right"
  ..- attr(*, "complete")= logi FALSE
  ..- attr(*, "validate")= logi TRUE
 $ coordinates:Classes 'CoordCartesian', 'Coord', 'ggproto', 'gg' <ggproto object: Class CoordCartesian, Coord, gg>
    aspect: function
    backtransform_range: function
    clip: on
    default: FALSE
    distance: function
    expand: TRUE
    is_free: function
    is_linear: function
    labels: function
    limits: list
    modify_scales: function
    range: function
    render_axis_h: function
    render_axis_v: function
    render_bg: function
    render_fg: function
    setup_data: function
    setup_layout: function
    setup_panel_guides: function
    setup_panel_params: function
    setup_params: function
    train_panel_guides: function
    transform: function
    super:  <ggproto object: Class CoordCartesian, Coord, gg> 
 $ facet      :Classes 'FacetNull', 'Facet', 'ggproto', 'gg' <ggproto object: Class FacetNull, Facet, gg>
    compute_layout: function
    draw_back: function
    draw_front: function
    draw_labels: function
    draw_panels: function
    finish_data: function
    init_scales: function
    map_data: function
    params: list
    setup_data: function
    setup_params: function
    shrink: TRUE
    train_scales: function
    vars: function
    super:  <ggproto object: Class FacetNull, Facet, gg> 
 $ plot_env   :<environment: 0x0000025b30e02ea0> 
 $ layout     :Classes 'Layout', 'ggproto', 'gg' <ggproto object: Class Layout, gg>
    coord: NULL
    coord_params: list
    facet: NULL
    facet_params: list
    finish_data: function
    get_scales: function
    layout: NULL
    map_position: function
    panel_params: NULL
    panel_scales_x: NULL
    panel_scales_y: NULL
    render: function
    render_labels: function
    reset_scales: function
    resolve_label: function
    setup: function
    setup_panel_guides: function
    setup_panel_params: function
    train_position: function
    super:  <ggproto object: Class Layout, gg> 
 $ labels     :List of 8
  ..$ title   : chr "Histogram and theoretical densities"
  ..$ y       : chr "Density"
  ..$ x       : chr "data"
  ..$ group   : chr "ind"
  ..$ colour  : chr "ind"
  ..$ linetype: chr "ind"
  ..$ size    : chr "ind"
  ..$ weight  : chr "weight"
  .. ..- attr(*, "fallback")= logi TRUE
 - attr(*, "class")= chr [1:2] "gg" "ggplot"
Warning message:
The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
ℹ The deprecated feature was likely used in the fitdistrplus package.
  Please report the issue at
  <https://github.com/lbbe-software/fitdistrplus/issues>. 
> 
> #test ylim argument
> denscomp(list(fitW, fitln, fitg), probability=TRUE, ylim=c(0, .05))
> denscomp(list(fitW, fitln, fitg), probability=FALSE, ylim=c(0, 100))
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(fitW, fitln, fitg), probability=TRUE, ylim=c(0, .05), plotstyle = "ggplot")
+   denscomp(list(fitW, fitln, fitg), probability=FALSE, ylim=c(0, 100), plotstyle = "ggplot")
+ }
> 
> #test xlim, legend, main, demp
> denscomp(list(fitW, fitln, fitg), legendtext=c("Weibull","lognormal","gamma"),
+          main="ground beef fits",xlab="serving sizes (g)",
+          ylab="F",xlim = c(0,250), xlegend = "topright", demp=TRUE)
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(fitW, fitln, fitg), legendtext=c("Weibull","lognormal","gamma"),
+            main="ground beef fits",xlab="serving sizes (g)",
+            ylab="F",xlim = c(0,250), xlegend = "topright", demp=TRUE, plotstyle = "ggplot")
+ }
> 
> 
> # (2) Plot lognormal distributions fitted by
> # maximum goodness-of-fit estimation
> # using various distances (data plotted in log scale)
> #
> data(endosulfan)
> ATV <-subset(endosulfan, group == "NonArthroInvert")$ATV
> flnMGEKS <- fitdist(ATV,"lnorm",method="mge",gof="KS")
> flnMGEAD <- fitdist(ATV,"lnorm",method="mge",gof="AD")
> flnMGEADL <- fitdist(ATV,"lnorm",method="mge",gof="ADL")
> flnMGEAD2L <- fitdist(ATV,"lnorm",method="mge",gof="AD2L")
> llfit <- list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L)
> 
> denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+          main="fits of a lognormal dist. using various GOF dist.")
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+            main="fits of a lognormal dist. using various GOF dist.", plotstyle = "ggplot")
+ }
> 
> denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+          main="fits of a lognormal dist. using various GOF dist.",
+          legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"))
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+            main="fits of a lognormal dist. using various GOF dist.",
+            legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"), plotstyle = "ggplot")
+ }
> 
> denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+          main="fits of a lognormal dist. using various GOF dist.",
+          legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"),
+          fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"))
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
+            main="fits of a lognormal dist. using various GOF dist.",
+            legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"),
+            fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"), plotstyle = "ggplot")
+ }
> 
> denscomp(llfit,
+          main="fits of a lognormal dist. using various GOF dist.",
+          legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"),
+          fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"),
+          datacol="grey")
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(llfit,
+            main="fits of a lognormal dist. using various GOF dist.",
+            legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"),
+            fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"),
+            datacol="grey", plotstyle = "ggplot")
+ }
> 
> denscomp(flnMGEKS, xlim=c(10,100000))
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize)
+   denscomp(flnMGEKS, xlim=c(10,100000), plotstyle = "ggplot")
> 
> 
> # (3)
> #
> #
> 
> if (visualize)
+ {
+   x1 <- c(6.4,13.3,4.1,1.3,14.1,10.6,9.9,9.6,15.3,22.1,13.4,13.2,8.4,6.3,8.9,5.2,10.9,14.4)
+   x <- seq(0, 1.1*max(x1), length=100)
+   
+   dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
+   pgumbel <- function(x,a,b) exp(-exp((a-x)/b))
+   
+   f1 <- mledist(x1,"norm")
+   f2 <- mledist(x1,"gumbel", start = list(a = 10, b = 5))
+   f3 <- mledist(x1, "exp")
+   
+   # graph 1
+   hist(x1, 10, prob=TRUE)
+   lines(x, dnorm(x, f1$estimate[1], f1$estimate[2]), col="red")
+   lines(x, dgumbel(x, f2$estimate[1], f2$estimate[2]), col="green")
+   lines(x, dexp(x, f3$estimate[1]), col="blue")
+   legend("topright", lty=1, leg = c("Normal", "Gumbel", "Exp"), col = c("red", "green", "blue"))
+   
+   # graph 2
+   f1 <- fitdist(x1,"norm")
+   f2 <- fitdist(x1,"gumbel", start = list(a = 10, b = 5))
+   f3 <- fitdist(x1, "exp")
+   denscomp(list(f1, f2, f3), xlim = c(0, 30), fitlty = 1, legendtext = c("Normal","Gumbel","Exp"))
+   
+   # graph 3
+   if (requireNamespace ("ggplot2", quietly = TRUE))
+     denscomp(list(f1, f2, f3), xlim = c(0, 30), fitlty = 1, legendtext = c("Normal","Gumbel","Exp"), breaks = 12, plotstyle = "ggplot")
+   
+   
+ }
> 
> # (4) normal mixture
> #
> if (visualize)
+ {
+   
+   #mixture of two normal distributions
+   #density
+   dnorm2 <- function(x, poid, m1, s1, m2, s2)
+     poid*dnorm(x, m1, s1) + (1-poid)*dnorm(x, m2, s2)
+   #numerical approximate quantile function
+   qnorm2 <- function(p, poid, m1, s1, m2, s2)
+   {
+     L2 <- function(x, prob)
+       (prob - pnorm2(x, poid, m1, s1, m2, s2))^2
+     sapply(p, function(pr) optimize(L2, c(-1000, 1000), prob=pr)$minimum)
+   }
+   #distribution function
+   pnorm2 <- function(q, poid, m1, s1, m2, s2)
+     poid*pnorm(q, m1, s1) + (1-poid)*pnorm(q, m2, s2)
+   
+   
+   #basic normal distribution
+   set.seed(1234)
+   x2 <- c(rnorm(nsample, 5),  rnorm(nsample, 10))
+   #MLE fit
+   fit1 <- fitdist(x2, "norm2", "mle", start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2),
+                   lower=c(0, 0, 0, 0, 0))
+   fit2 <- fitdist(x2, "norm2", "qme", probs=c(1/6, 1/4, 1/3, 1/2, 2/3),
+                   start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2),
+                   lower=c(0, 0, 0, 0, 0), upper=c(1/2, Inf, Inf, Inf, Inf))
+   fit3 <- fitdist(x2, "norm2", "mge", gof="AD",
+                   start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2),
+                   lower=c(0, 0, 0, 0, 0), upper=c(1/2, Inf, Inf, Inf, Inf))
+   
+   denscomp(list(fit1, fit2, fit3))
+   if (requireNamespace ("ggplot2", quietly = TRUE) & visualize)
+     denscomp(list(fit1, fit2, fit3), plotstyle = "ggplot")
+ }
> 
> # (5) large data
> #
> if (visualize)
+ {
+   n <- 1e4
+   x <- rnorm(n)
+   f <- fitdist(x, "norm")
+   
+   denscomp(f)
+   denscomp(f, demp=TRUE)
+   if (requireNamespace ("ggplot2", quietly = TRUE)) {
+     denscomp(f, plotstyle = "ggplot")
+     denscomp(f, demp=TRUE, plotstyle = "ggplot")
+   }
+ }
> 
> 
> # (6) graphical parameters
> #
> if (visualize)
+ {
+   # 'graphics' plot style
+   denscomp(list(fit1, fit2, fit3), plotstyle = "gr")
+   denscomp(list(fit1, fit2, fit3), title = "Fitted distribution")
+   denscomp(list(fit1, fit2, fit3), main = "Fitted distribution", addlegend = F, demp = T, dempcol = "purple")
+   
+   # 'ggplot' plot style
+   if (requireNamespace ("ggplot2", quietly = TRUE)) {
+     denscomp(list(fit1, fit2, fit3), plotstyle = "gg")
+     denscomp(list(fit1, fit2, fit3), plotstyle = "ggplot", breaks = 20, pro = F)
+     dcomp <- denscomp(list(fit1, fit2, fit3), plotstyle = "gg", demp = T)
+     dcomp + ggplot2::theme_minimal() + ggplot2::ggtitle("Histogram and\ntheoretical densities")
+     dcomp + ggplot2::guides(colour = ggplot2::guide_legend("Fitted distribution"), linetype = ggplot2::guide_legend("Fitted distribution")) 
+   }
+ }
> 
> 
> 
> # (7) test legend labels
> #
> if (visualize)
+ {
+   serving <- groundbeef$serving
+   fitW <- fitdist(serving,"weibull")
+   fitW2 <- fitdist(serving,"weibull", method="qme", probs=c(1/3,2/3))
+   fitW3 <- fitdist(serving,"weibull", method="qme", probs=c(1/2,2/3))
+   fitln <- fitdist(serving,"lnorm")
+   fitg <- fitdist(serving,"gamma")
+   
+   denscomp(list(fitW, fitln, fitg)) #distrib
+   denscomp(list(fitW, fitW2, fitln, fitg)) #distrib+method
+   denscomp(list(fitW, fitW2, fitW3, fitln, fitg)) #distrib+method+num
+   if (requireNamespace ("ggplot2", quietly = TRUE))
+     denscomp(list(fitW, fitW2, fitW3, fitln, fitg), plotstyle = "ggplot") #distrib+method+num
+ }
> 
> # (8) discrete distrib
> #
> 
> x <- c(rpois(nsample, 5), rbinom(nsample, 12, 2/3))
> fpois <- fitdist(x, "pois")
> fgeo <- fitdist(x, "geom")
> fnbinom <- fitdist(x, "nbinom")
> 
> par(mar=c(4,4,2,1))
> denscomp(list(fpois, fnbinom, fgeo), probability = TRUE)
> denscomp(list(fpois, fnbinom, fgeo), probability = FALSE)
> denscomp(list(fpois, fnbinom, fgeo), fittype="o")
> denscomp(list(fpois, fnbinom, fgeo), fittype="p")
> # 'ggplot' plot style
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
+   denscomp(list(fpois, fnbinom, fgeo), plotstyle="ggplot", probability = TRUE)
+   denscomp(list(fpois, fnbinom, fgeo), plotstyle="ggplot", probability = FALSE)
+   denscomp(list(fpois, fnbinom, fgeo), fittype="o", plotstyle="ggplot")
+   denscomp(list(fpois, fnbinom, fgeo), fittype="p", plotstyle="ggplot")
+ }
> 
> # test the call to any()
> fpois$discrete <- fnbinom$discrete <- FALSE
> denscomp(list(fpois, fnbinom, fgeo))
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize)
+   denscomp(list(fpois, fnbinom, fgeo), plotstyle="ggplot")
> #test the forced usage
> fgeo$discrete <- FALSE
> denscomp(list(fpois, fnbinom, fgeo), discrete=TRUE)
> if (requireNamespace ("ggplot2", quietly = TRUE) & visualize)
+   denscomp(list(fpois, fnbinom, fgeo), discrete=TRUE, plotstyle="ggplot")
> 
> if (visualize)
+ {
+   x <- c(rpois(nsample, 30), rbinom(nsample, 12, 2/3))
+   fpois <- fitdist(x, "pois")
+   fgeo <- fitdist(x, "geom")
+   fnbinom <- fitdist(x, "nbinom")
+   
+   #3 types of plot of probability mass function
+   par(mar=c(4,4,2,1))
+   denscomp(list(fpois, fnbinom, fgeo)) #fittype = "l"
+   denscomp(list(fpois, fnbinom, fgeo), fittype = "p")
+   denscomp(list(fpois, fnbinom, fgeo), fittype = "o")
+   if (requireNamespace ("ggplot2", quietly = TRUE)) {
+     denscomp(list(fpois, fnbinom, fgeo), plotstyle="ggplot") #fittype = "l"
+     denscomp(list(fpois, fnbinom, fgeo), fittype = "p", plotstyle="ggplot")
+     denscomp(list(fpois, fnbinom, fgeo), fittype = "o", plotstyle="ggplot")
+   }
+ }
> 
> # (9) examples with user specified regular of irregular breaks in the histogram
> #     in probability or not
> #
> if (visualize)
+ {
+   # two plots with user specified regular breaks in probability or not
+   # hist(serving, breaks = seq(0,200,50))
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            breaks = seq(0,200,50))
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            probability = FALSE, breaks = seq(0,200,50))
+   # with ggplot2
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            plotstyle = "ggplot", breaks = seq(0,200,50))
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            probability = FALSE, plotstyle = "ggplot", breaks = seq(0,200,50))
+   
+   
+   # two plots with irregular breaks in probability or not
+   # hist(serving, breaks = c(0, 20, 50, 100, 200, 300))
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            breaks = c(0, 20, 50, 100, 200, 300))
+   # hist(serving, breaks = c(0, 20, 50, 100, 200, 300), probability = FALSE)
+   try(denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+                main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+                breaks = c(0, 20, 50, 100, 200, 300), probability = FALSE))
+   # with ggplot2
+   denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+            main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+            breaks = c(0, 20, 50, 100, 200, 300), plotstyle = "ggplot")
+   ##### ggplot2 does not take into account non-equidistant breaks !!!!!!!!!!!!!!!!
+   try(denscomp(list(fitW, fitln, fitg), addlegend = FALSE,
+                main = "ground beef fits", xlab = "serving sizes (g)", xlim = c(0, 250), 
+                breaks = c(0, 20, 50, 100, 200, 300), 
+                probability = FALSE, plotstyle = "ggplot"))
+ }
> 
> 
> # (10) fitlty, fitlwd for discrete
> x <- c(rpois(nsample, 30), rbinom(nsample, 12, 2/3))
> fpois <- fitdist(x, "pois")
> fgeo <- fitdist(x, "geom")
> fnbinom <- fitdist(x, "nbinom")
> denscomp(list(fpois, fnbinom, fgeo), fitlty = 2, fitlwd = 3:1)
> denscomp(list(fpois, fnbinom, fgeo), fittype = "o", fitlwd = 3:1)
> if (requireNamespace ("ggplot2", quietly = TRUE)) {
+   denscomp(list(fpois, fnbinom, fgeo), plotstyle = "ggplot", fitlty = 2, fitlwd = 3:1)
+   denscomp(list(fpois, fnbinom, fgeo), plotstyle = "ggplot", fittype = "o", fitlty = 1, fitlwd = 3:1)
+ }
> 
> # (11) fitlty, fitlwd for non discrete
> denscomp(list(fitW, fitln, fitg), fitlty = 1, fitlwd = 3:1)
> denscomp(list(fitW, fitln, fitg), fitlty = 1, fitlwd = 1:3, fitcol = c(1:2, 7))
> if (requireNamespace ("ggplot2", quietly = TRUE)) {
+   denscomp(list(fitW, fitln, fitg), plotstyle = "ggplot", fitlty = 1, fitlwd = 3:1)
+   denscomp(list(fitW, fitln, fitg), plotstyle = "ggplot", fitlty = 1, fitlwd = 1:3, fitcol = c(1:2, 7))
+ }
> 
> proc.time()
   user  system elapsed 
   3.96    0.51    4.48