R Under development (unstable) (2025-01-06 r87534 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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 <- 1e3
> nsample <- 10
> 
> #### (1) Gamma example ####
> 
> truetheta <- c("alpha"=0.19, "beta"=5.18)
> x <- rgamma(nsample, truetheta["alpha"], truetheta["beta"])
> f1 <- mledist(x, "gamma", calcvcov = TRUE)
> f1$vcov
            shape       rate
shape 0.008117369  0.1676069
rate  0.167606943 14.4665363
> f1$estimate
    shape      rate 
0.2581476 5.3302163 
> infoFisher <- function(alpha, beta)
+ {
+   cbind(c(trigamma(alpha), -1/beta),
+         c(-1/beta, alpha/beta))
+ }
> solve(infoFisher(0.19, 5.18))/nsample
           [,1]       [,2]
[1,] 0.00357563 0.01881911
[2,] 0.01881911 2.82536371
> 
> 
> if(FALSE)
+ {
+   #check with MASS::fitdistr()
+   
+   mledist(rgamma(1e2, truetheta["alpha"], truetheta["beta"]), "gamma", calcvcov = TRUE)$vcov
+   mledist(rgamma(1e3, truetheta["alpha"], truetheta["beta"]), "gamma", calcvcov = TRUE)$vcov
+   mledist(rgamma(1e4, truetheta["alpha"], truetheta["beta"]), "gamma", calcvcov = TRUE)$vcov
+   
+   
+   MASS::fitdistr(rgamma(1e2, truetheta["alpha"], truetheta["beta"]), "gamma")$vcov
+   MASS::fitdistr(rgamma(1e3, truetheta["alpha"], truetheta["beta"]), "gamma")$vcov
+   MASS::fitdistr(rgamma(1e4, truetheta["alpha"], truetheta["beta"]), "gamma")$vcov
+   
+   
+ }
> 
> 
> # (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.3334977  2.7598598 22.3043053 

$convergence
[1] 0

$value
[1] 0.1049306

$hessian
          rate low upp
rate 0.4495647   0   0
low  0.0000000   0   0
upp  0.0000000   0   0

$optim.function
[1] "optim"

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

$fix.arg
NULL

$fix.arg.fun
NULL

$weights
NULL

$counts
function gradient 
      82       NA 

$optim.message
NULL

$loglik
[1] -10.49306

$vcov
NULL

> 
> proc.time()
   user  system elapsed 
   1.17    0.17    1.32