correctly normalize the asymptotic vcov

lbbe-software · Dec 5, 2024 · 6a5fa41 · 6a5fa41
1 parent f23b150
commit 6a5fa41
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 9 deletions.
diff --git a/R/mledist.R b/R/mledist.R
@@ -404,7 +404,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
     if(calcvcov && !is.null(opt$hessian))
     {
       #see R/util-mledist-vcov.R
-      varcovar <- mle.vcov(opt$hessian)
+      varcovar <- mle.vcov(opt$hessian, data.size)
       #add names
       if(!is.null(varcovar))
         colnames(varcovar) <- rownames(varcovar) <- names(opt$par)
@@ -445,7 +445,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
       names(opt$par) <- names(vstart)
     if(calcvcov && !is.null(opt$hessian))
     {
-      varcovar <- mle.vcov(opt$hessian)
+      varcovar <- mle.vcov(opt$hessian, data.size)
       #add names
       if(!is.null(varcovar))
         colnames(varcovar) <- rownames(varcovar) <- names(opt$par)

diff --git a/R/util-mledist-vcov.R b/R/util-mledist-vcov.R
@@ -23,12 +23,12 @@
 ### 
 
 
-#compute H^{-1}
-mle.vcov <- function(myhessian)
+#compute H^{-1}/n
+mle.vcov <- function(myhessian, nsample)
 {
   if(all(!is.na(myhessian)) && qr(myhessian)$rank == NCOL(myhessian))
   {
-    res <- solve(myhessian)
+    res <- solve(myhessian)/nsample
   }else
   {
     res <- NULL

diff --git a/R/util-mmedist-vcov.R b/R/util-mmedist-vcov.R
@@ -24,7 +24,7 @@
 ### 
 
 
-#compute J^{-1} A J^{-T}, see below
+#compute J^{-1} A J^{-T}/n, see below
 mme.vcov <- function(par, fix.arg, order, obs, mdistnam, memp, weights, 
                      epsilon = sqrt(.Machine$double.eps), echo=FALSE)
 {
@@ -39,7 +39,7 @@ mme.vcov <- function(par, fix.arg, order, obs, mdistnam, memp, weights,
   {  
     Jinv <- solve(Jmat)
     Amat <- mme.Ahat(par, fix.arg, order, obs, mdistnam, memp, weights)
-    res <- Jinv %*% Amat %*% t(Jinv)
+    res <- Jinv %*% Amat %*% t(Jinv) / length(obs)
   }else 
     res <- NULL
   if(echo)

diff --git a/tests/t-mledist-asymptotic-vcov.R b/tests/t-mledist-asymptotic-vcov.R
@@ -1,13 +1,38 @@
 require("fitdistrplus")
-nsample <- 1e6
+nsample <- 1e3
 nsample <- 10
 
 #### (1) Gamma example ####
 
-truetheta <- c("alpha"=3, "beta"=1/2)
+truetheta <- c("alpha"=0.19, "beta"=5.18)
 x <- rgamma(nsample, truetheta["alpha"], truetheta["beta"])
 f1 <- mledist(x, "gamma", calcvcov = TRUE)
 f1$vcov
+f1$estimate
+infoFisher <- function(alpha, beta)
+{
+  cbind(c(trigamma(alpha), -1/beta),
+        c(-1/beta, alpha/beta))
+}
+solve(infoFisher(0.19, 5.18))/nsample
+
+
+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
 #

diff --git a/tests/t-mmedist-asymptotic-vcov.R b/tests/t-mmedist-asymptotic-vcov.R
@@ -15,6 +15,13 @@ if(FALSE)
   fitdistrplus:::mme.vcov(as.numeric(truetheta), fix.arg=NULL, order=1:2, obs=x, mdistnam=mgamma, memp, weights=NULL)
 }
 
+if(FALSE)
+{
+  mmedist(rgamma(1e2, truetheta["alpha"], truetheta["beta"]), "gamma", order=1:2, calcvcov = TRUE)$vcov
+  mmedist(rgamma(1e3, truetheta["alpha"], truetheta["beta"]), "gamma", order=1:2, calcvcov = TRUE)$vcov
+  mmedist(rgamma(1e4, truetheta["alpha"], truetheta["beta"]), "gamma", order=1:2, calcvcov = TRUE)$vcov
+}
+
 
 # (2) fit a Pareto distribution
 #