The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
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[@]dogado.de.

Bootstrap implementation

Jaime Mosquera

2022-12-09

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:     43112.561:  1.00000
#>   1:     32427.035:  2.00000
#>   2:     32176.585:  1.91483
#>   3:     32027.228:  1.75021
#>   4:     32016.630:  1.78713
#>   5:     32016.434:  1.78285
#>   6:     32016.434:  1.78270
#>   7:     32016.434:  1.78270
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Hessian from optim 
#> _______________________________________________________________
#>        AIC      BIC
#>   64032.87 64032.87
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   5.94592    0.04204   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,] 565.708

## Standard errors
print(theta_1$fit$StdE)
#> [1] 0.04204398
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:     43112.561:  1.00000
#>   1:     32427.035:  2.00000
#>   2:     32176.585:  1.91483
#>   3:     32027.228:  1.75021
#>   4:     32016.630:  1.78713
#>   5:     32016.434:  1.78285
#>   6:     32016.434:  1.78270
#>   7:     32016.434:  1.78270
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
#>   64032.87 64032.87
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   5.94592    0.04582   129.8   <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,] 565.708

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

Notice that Standard Errors calculated with optim (\(0.042044\)) and those calculated with bootstrap implementation (\(0.045823\)) 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.
Health stats visible at Monitor.