Bootstrap implementation
Jaime Mosquera
2025-09-05
Source:vignettes/Bootstrap_implementation.Rmd
Bootstrap_implementation.Rmd
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:
where is the observed Fisher Information Matrix. Hence, the standard errors can be calculated as the square root of the diagonal elements of matrix , as follows:
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:
- If
StdE_Method = optim
, it is estimated through theoptim
function (with optionhessian = TRUE
under the hood inmaxlogL
ormaxlogLreg
function). - If the previous implementation fails or if the user chooses
StdE_Method = numDeriv
, it is calculated withhessian
function from numDeriv package. - If both of the previous methods fail, then standard errors are
computed by bootstrapping with the function
bootstrap_maxlogL
.
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: 42998.545: 1.00000
#> 1: 32411.605: 2.00000
#> 2: 32155.569: 1.91367
#> 3: 32003.818: 1.74712
#> 4: 31992.756: 1.78488
#> 5: 31992.548: 1.78047
#> 6: 31992.548: 1.78031
#> 7: 31992.548: 1.78032
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> 63985.1 63985.1
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> sd 5.93173 0.04194 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,] 568.417
## Standard errors
print(theta_1$fit$StdE)
#> [1] 0.04194367
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: 42998.545: 1.00000
#> 1: 32411.605: 2.00000
#> 2: 32155.569: 1.91367
#> 3: 32003.818: 1.74712
#> 4: 31992.756: 1.78488
#> 5: 31992.548: 1.78047
#> 6: 31992.548: 1.78031
#> 7: 31992.548: 1.78032
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
#> 63985.1 63985.1
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> sd 5.9317 0.0401 147.9 <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,] 568.417
## Standard errors
print(theta_2$fit$StdE)
#> [1] 0.04010111
print(theta_2$outputs$StdE_Method)
#> [1] "Bootstrap"
Notice that Standard Errors calculated with optim
()
and those calculated with bootstrap implementation
()
are approximately equals, but no identical.