The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

Bootstrap implementation

Jaime Mosquera

2025-09-08

Computation of standard errors

Objects of maxlogL class (outputs from maxlogL and maxlogLreg) stores the estimated parameters of probability density/mass functions by Maximum Likelihood. The variance-covariance matrix is computed from Fisher information matrix, which is obtained by means of the Inverse Hessian matrix of estimators:

\[\begin{equation} Var(\hat{\boldsymbol{\theta}}) = \mathcal{J}^{-1}(\hat{\boldsymbol{\theta}}) = C(\hat{\boldsymbol{\theta}}), \end{equation}\]

where \(\mathcal{J}(\hat{\boldsymbol{\theta}})\) is the observed Fisher Information Matrix. Hence, the standard errors can be calculated as the square root of the diagonal elements of matrix \(C\), as follows:

\[\begin{equation} SE(\hat{\boldsymbol{\theta}}) = \sqrt{C_{jj}(\hat{\boldsymbol{\theta}})}, \end{equation}\]

To install the package, type the following commands:

if (!require('devtools')) install.packages('devtools')
devtools::install_github('Jaimemosg/EstimationTools', force = TRUE)

In EstimationTools Hessian matrix is computed in the following way:

Additionally, EstimationTools allows implementation of bootstrap for standard error estimation, even if the Hessian computation does not fail.

Standard Error with maxlogL function

Lets fit the following distribution:

\[ \begin{aligned} X &\sim N(\mu, \:\sigma^2) \\ \mu &= 160 \quad (\verb|mean|) \\ \sigma &= 6 \quad (\verb|sd|) \end{aligned} \]

The following chunk illustrates the fitting with Hessian computation via optim:

library(EstimationTools)

x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                   link = list(over = "sd", fun = "log_link"),
                   fixed = list(mean = 160))
#>   0:     43580.164:  1.00000
#>   1:     32490.318:  2.00000
#>   2:     32261.854:  1.91946
#>   3:     32122.237:  1.76266
#>   4:     32113.373:  1.79628
#>   5:     32113.222:  1.79250
#>   6:     32113.221:  1.79238
#>   7:     32113.221:  1.79238
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Hessian from optim 
#> _______________________________________________________________
#>        AIC      BIC
#>   64226.44 64226.44
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   6.00375    0.04245   141.4   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators 
#> ---

## Hessian
print(theta_1$fit$hessian)
#>          [,1]
#> [1,] 554.8626

## Standard errors
print(theta_1$fit$StdE)
#> [1] 0.04245289
print(theta_1$outputs$StdE_Method)
#> [1] "Hessian from optim"

Note that Hessian was computed with no issues. Now, lets check the aforementioned feature in maxlogL: the user can implement bootstrap algorithm available in bootstrap_maxlogL function. To illustrate this, we are going to create another object theta_2:

# Bootstrap
theta_2 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                   link = list(over = "sd", fun = "log_link"),
                   fixed = list(mean = 160))
#>   0:     43580.164:  1.00000
#>   1:     32490.318:  2.00000
#>   2:     32261.854:  1.91946
#>   3:     32122.237:  1.76266
#>   4:     32113.373:  1.79628
#>   5:     32113.222:  1.79250
#>   6:     32113.221:  1.79238
#>   7:     32113.221:  1.79238
bootstrap_maxlogL(theta_2, R = 200)
#> 
#> ...Bootstrap computation of Standard Error. Please, wait a few minutes...
#> 
#>  --> Done <---
summary(theta_2)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Bootstrap 
#> _______________________________________________________________
#>        AIC      BIC
#>   64226.44 64226.44
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   6.00375    0.04256   141.1   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators 
#> ---

## Hessian
print(theta_2$fit$hessian)
#>          [,1]
#> [1,] 554.8626

## Standard errors
print(theta_2$fit$StdE)
#> [1] 0.04255708
print(theta_2$outputs$StdE_Method)
#> [1] "Bootstrap"

Notice that Standard Errors calculated with optim (\(0.042453\)) and those calculated with bootstrap implementation (\(0.042557\)) are approximately equals, but no identical.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.