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 > nsample <- 10 > > # (1) basic fit of a normal distribution with moment matching estimation > # > > set.seed(1234) > x1 <- rnorm(n=nsample) > mmedist(x1,"norm") $estimate mean sd -0.3831574 0.9446870 $convergence [1] 0 $value NULL $hessian NULL $optim.function NULL $opt.meth NULL $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts NULL $optim.message NULL $loglik [1] -13.62037 $method [1] "closed formula" $order [1] 1 2 $memp NULL $vcov NULL > > try(mmedist(x1,"norm", fix.arg=list(mean=0))) Error in mmedist(x1, "norm", fix.arg = list(mean = 0)) : argument fix.arg cannot be used when a closed formula is used. > > # (2) fit a discrete distribution (Poisson) > # > > set.seed(1234) > x2 <- rpois(n=nsample, lambda = 2) > mmedist(x2,"pois") $estimate lambda 1.7 $convergence [1] 0 $value NULL $hessian NULL $optim.function NULL $opt.meth NULL $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts NULL $optim.message NULL $loglik [1] -15.31626 $method [1] "closed formula" $order [1] 1 $memp NULL $vcov NULL > > # (3) fit a finite-support distribution (beta) > # > > set.seed(1234) > x3 <- rbeta(n=nsample, shape1=5, shape2=10) > mmedist(x3,"beta") $estimate shape1 shape2 3.956081 9.512364 $convergence [1] 0 $value NULL $hessian NULL $optim.function NULL $opt.meth NULL $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts NULL $optim.message NULL $loglik [1] 6.939847 $method [1] "closed formula" $order [1] 1 2 $memp NULL $vcov NULL > > > > # (4) 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 + mmedist(x4, "pareto", order=c(1, 2), memp=memp, start=c(shape=10, scale=10), + lower=1, upper=Inf) + mmedist(x4, "pareto", order=1, memp=memp, start=list(shape=10), fix.arg=list(scale=1.5), + lower=2, upper=Inf) + mmedist(x4, "pareto", order=1, memp=memp, start=function(x) list(shape=10), fix.arg=list(scale=1.5), + lower=2, upper=Inf) + mmedist(x4, "pareto", order=1, memp=memp, start=list(shape=10), fix.arg=function(x) list(scale=1.5), + lower=2, upper=Inf) + + #weights + memp2 <- function(x, order, weights) sum(x^order * weights)/sum(weights) + w <- rep(1, length(x4)) + w[x4 < 1] <- 2 + mmedist(x4, "pareto", order=c(1, 2), memp=memp2, weights=w, + start=list(shape=10, scale=10), lower=1, upper=Inf) + + #fit + data(danishuni) + fparedanishMME <- mmedist(danishuni$Loss, "pareto", order=1:2, + memp=memp, start=c(shape=10, scale=10), lower=2+1e-6, upper=Inf) + c(theo = mpareto(1, fparedanishMME$estimate[1], fparedanishMME$estimate[2]), + emp = memp(danishuni$Loss, 1)) + c(theo = mpareto(2, fparedanishMME$estimate[1], fparedanishMME$estimate[2]), + emp = memp(danishuni$Loss, 2)) + + } 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 theo emp 83.80216 83.80216 Warning messages: 1: In mmedist(x4, "pareto", order = 1, memp = memp, start = list(shape = 10), : The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B. 2: In mmedist(x4, "pareto", order = 1, memp = memp, start = function(x) list(shape = 10), : The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B. 3: In mmedist(x4, "pareto", order = 1, memp = memp, start = list(shape = 10), : The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B. 4: In mmedist(x4, "pareto", order = c(1, 2), memp = memp2, weights = w, : weights are not taken into account in the default initial values > > > # (5) fit a lognormal distribution > # > x3 <- rlnorm(n=nsample) > f1 <- mledist(x3, "lnorm") #previously mmedist was the same as mledist > f2 <- mmedist(x3, "lnorm") > n <- length(x3) > s2 <- log(1+var(x3)/mean(x3)^2*(n-1)/n) > mu <- log(mean(x3)) - s2/2 > cbind(c(mu, s2), f2$estimate) [,1] [,2] meanlog -0.1123774 -0.1123774 sdlog 1.4016705 1.1839216 > > > c(truestim=exp(mu+s2/2), + jensen=as.numeric(exp(f1$estimate["meanlog"]+f1$estimate["sdlog"]^2/2)), + emp=mean(x3)) truestim jensen emp 1.801209 1.399251 1.801209 > > c(truestim=exp(2*mu+s2)*(exp(s2)-1), + jensen=as.numeric(exp(2*f1$estimate["meanlog"]+f1$estimate["sdlog"]^2)*(exp(f1$estimate["sdlog"]^2)-1)), + emp=var(x3)*(n-1)/n) truestim jensen emp 9.934141 2.710276 9.934141 > > > # (6) test error messages > # > > mnorm3 <- dnorm3 <- function(x, a) + "NA" > x <- rnorm(nsample) > > #should get a one-line error > res <- mmedist(x, "norm3", start=list(a=1), order=1, memp=function(x, order) mean(x)) > #as in > attr(try(log("a"), silent=TRUE), "condition") > > > try(mmedist(c(serving, "a"), "gamma")) Error : object 'serving' not found > try(mmedist(c(serving, NA), "gamma")) Error : object 'serving' not found > try(mmedist(c(serving, Inf), "gamma")) Error : object 'serving' not found > try(mmedist(c(serving, -Inf), "gamma")) Error : object 'serving' not found > try(mmedist(c(serving, NaN), "gamma")) Error : object 'serving' not found > > > # (7) fit of a normal distribution with weighted moment matching estimation > # > > n <- nsample > w <- c(rep(1, n/2), rep(10, n/2)) > mmedist(x1, "norm", weights=w)$estimate mean sd -0.4083605 0.5960473 Warning message: In mmedist(x1, "norm", weights = w) : weights are not taken into account in the default initial values > > #check > sum(w*x1)/sum(w) [1] -0.4083605 > fitdistrplus:::wtdmean(x1, w) [1] -0.4083605 > sum(w*(x1-sum(w*x1)/sum(w))^2)/sum(w) [1] 0.3488129 > fitdistrplus:::wtdvar(x1, w) [1] 0.3552724 > > > mmedist(exp(x1), "lnorm", weights=w)$estimate meanlog sdlog -0.4199772 0.6276551 Warning message: In mmedist(exp(x1), "lnorm", weights = w) : weights are not taken into account in the default initial values > > #test non integer weights > try(mmedist(x1, "norm", weights=rep(1/3, length(x1)))) Error in mmedist(x1, "norm", weights = rep(1/3, length(x1))) : weights should be a vector of (strictly) positive integers > try(mmedist(1:10, "pois", weights=c(rep(1, 9), 1.001), start=list(lambda=mean(x)))) Error in mmedist(1:10, "pois", weights = c(rep(1, 9), 1.001), start = list(lambda = mean(x))) : weights should be a vector of (strictly) positive integers > try(mmedist(1:10, "pois", weights=c(rep(1, 9), 1.0000001), start=list(lambda=mean(x)))) Error in mmedist(1:10, "pois", weights = c(rep(1, 9), 1.0000001), start = list(lambda = mean(x))) : weights should be a vector of (strictly) positive integers > > > # (8) fit of a neg binom distribution with weighted moment matching estimation > # > > x4 <- rnbinom(nsample, 5, 1/2) > table(x4) x4 0 1 2 4 5 6 7 8 2 1 1 1 2 1 1 1 > w <- rep(1, length(x4)) > w[x4 > 5] <- 2 > > mmedist(x4, "nbinom", weights=w)$estimate size mu 5.284919 4.538462 Warning message: In mmedist(x4, "nbinom", weights = w) : weights are not taken into account in the default initial values > mmedist(x4, "nbinom", weights=NULL)$estimate size mu 3.840426 3.800000 > > # (9) relevant example for zero modified geometric distribution > # > > memp1 <- function(x, order) mean(x^order) > memp2 <- function(x, order, weights) sum(x^order * weights)/sum(weights) > > > rzmgeom <- function(n, p1, p2) + { + u <- rbinom(n, 1, 1-p1) #prob to get zero is p1 + u[u != 0] <- rgeom(sum(u != 0), p2)+1 + u + } > dzmgeom <- function(x, p1, p2) + { + p1 * (x == 0) + (1-p1)*dgeom(x-1, p2) + } > mcumulgeom <- function(order, prob) #E[N(N-1)..(N-o+1)] + { + (1-prob)^order * factorial(order) / prob^order + } > mgeom <- function(order, prob) + { + if(order == 0) + return(1) + m1 <- mcumulgeom(1, prob) + if(order == 1) + return(m1) + m2 <- mcumulgeom(2, prob) + m1 + if(order == 2) + return(m2) + m3 <- mcumulgeom(3, prob) + 3*m2 + 2*m1 + if(order == 3) + return(m3) + m4 <- mcumulgeom(4, prob) + 6*m3 + 11*m2 - 6*m1 + if(order == 4) + return(m4) + else + stop("not yet implemented") + } > if(FALSE) + { + for(j in 1:4) + print(c(memp1(rgeom(1e4, 1/6), j), mgeom(j, 1/6))) + } > > mzmgeom <- function(order, p1, p2) #raw moment + { + momj <- sapply(0:order, function(j) mgeom(j, p2)) + cj <- choose(order, 0:order) + (1-p1) * sum( cj * momj ) + } > if(FALSE) + { + for(j in 1:4) + print(c(memp1(rzmgeom(1e4, 1/3, 1/6), j), mzmgeom(j, 1/3, 1/6))) + + } > > > x5 <- rzmgeom(nsample, 1/3, 1/6) > w <- rep(1, length(x5)) > w[x5 > 5] <- 2 > > mmedist(x5, "zmgeom", order=1:2, memp=memp1, start=list(p1=mean(x5 == 0), p2=1/mean(x5[x5 > 0])), + lower=0.01, upper=0.99)$estimate p1 p2 0.3382355 0.2205882 > mmedist(x5, "zmgeom", order=1:2, memp=memp2, start=list(p1=mean(x5 == 0), p2=1/mean(x5[x5 > 0])), + weights=w)$estimate [1] NA NA Warning message: In mmedist(x5, "zmgeom", order = 1:2, memp = memp2, start = list(p1 = mean(x5 == : weights are not taken into account in the default initial values > > > mmedist(x5, "zmgeom", order=1:2, memp=memp1, start=list(p1=mean(x5 == 0), p2=1/mean(x5[x5 > 0])), + lower=0.01, upper=0.99)$loglik [1] -21.27197 > mmedist(x5, "zmgeom", order=1:2, memp=memp2, start=list(p1=mean(x5 == 0), p2=1/mean(x5[x5 > 0])), + weights=w)$loglik NULL Warning message: In mmedist(x5, "zmgeom", order = 1:2, memp = memp2, start = list(p1 = mean(x5 == : weights are not taken into account in the default initial values > > # (10) bounds > # > 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 + mmedist(x4, "pareto", order=c(1, 2), memp=memp, start=c(shape=10, scale=10), lower=1, upper=Inf, optim.method = "L-BFGS-B") #L-BFGS-B via optim + mmedist(x4, "pareto", order=c(1, 2), memp=memp, start=c(shape=10, scale=10), lower=1, upper=Inf, optim.method = "Nelder") #Nelder Mead via constrOptim + } $estimate shape scale 51.643906 8.808025 $convergence [1] 7 $value [1] 0.0001558761 $hessian NULL $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 4818 NA $optim.message [1] "Barrier algorithm ran out of iterations and did not converge" $loglik [1] 7.031066 $method [1] "Nelder" $order [1] 1 2 $memp function (x, order) mean(x^order) $vcov NULL > > proc.time() user system elapsed 4.23 0.56 4.76