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.

> require(fitdistrplus)
Loading required package: fitdistrplus
Loading required package: MASS
Loading required package: survival
> nsample <- 1e6
> nsample <- 10
> 
> #### (1) Gamma example ####
> 
> truetheta <- c("alpha"=3, "beta"=1/2)
> x <- rgamma(nsample, truetheta["alpha"], truetheta["beta"])
> f1 <- mledist(x, "gamma", calcvcov = TRUE)
> f1$vcov
          shape       rate
shape 0.7802361 0.11247195
rate  0.1124719 0.02061668
> 
> # (2) fit a Pareto distribution
> #
> 
> if(any(installed.packages()[, "Package"] == "actuar"))
+ {
+   require(actuar)
+   #simulate a sample
+   x4 <- rpareto(nsample, 6, 2)
+   
+   #empirical raw moment
+   memp <- function(x, order)
+     mean(x^order)
+   
+   #fit
+   res <- mledist(x4, "pareto", start=list(shape=10, scale=10), 
+           lower=1, upper=Inf, calcvcov = TRUE)
+ }
Loading required package: actuar

Attaching package: 'actuar'

The following objects are masked from 'package:stats':

    sd, var

The following object is masked from 'package:grDevices':

    cm

> 
> # (3) truncated distribution
> #
> 
> dtiexp <- function(x, rate, low, upp)
+ {
+   PU <- pexp(upp, rate=rate, lower.tail = FALSE)
+   PL <- pexp(low, rate=rate)
+   dexp(x, rate) * (x >= low) * (x <= upp) + PL * (x == low) + PU * (x == upp)
+ }
> ptiexp <- function(q, rate, low, upp)
+   pexp(q, rate) * (q >= low) * (q <= upp) + 1 * (q > upp)
> n <- 100; x <- pmax(pmin(rexp(n), 3), .5)
> 
> mledist(x, "tiexp", start=list(rate=3, low=0, upp=20), calcvcov = TRUE)
$estimate
      rate        low        upp 
 0.9326635 -0.5099761 23.0592451 

$convergence
[1] 0

$value
[1] 106.9494

$hessian
         rate low upp
rate 114.9611   0   0
low    0.0000   0   0
upp    0.0000   0   0

$optim.function
[1] "optim"

$optim.method
[1] "Nelder-Mead"

$fix.arg
NULL

$fix.arg.fun
NULL

$weights
NULL

$counts
function gradient 
      64       NA 

$optim.message
NULL

$loglik
[1] -106.9494

$vcov
NULL

> 
> proc.time()
   user  system elapsed 
   1.98    0.42    2.34