Integration and Differentiation

Install required packages

numDeriv

wants <- c("numDeriv")
has   <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])

Integration

Example

Half the normal density is to the left of the mean

integrate(dnorm, -Inf, 1, mean=1, sd=2)$value

[1] 0.5

pnorm(1, mean=1, sd=2)

[1] 0.5

Re-implement normal cumulative distribution function

Case study for implementing a cumulative distribution function which may not be available in closed form whereas the corresponding density function is.

Naive implementation of normal density - numerically unstable, only works within narrow parameter range.

dGauss <- function(x, mu=0, sigma=1) {
    (1/(sigma*sqrt(2*pi))) * exp(-0.5 * ((x-mu)/sigma)^2)
}

Implement normal cumulative distribution function by integrating density function.

pGauss <- function(x, mu=0, sigma=1) {
    integrate(dGauss, lower=-Inf, upper=x, mu=mu, sigma=sigma)$value
}

Check custom CDF against R implementation

pGauss(1, mu=0, sigma=1)

[1] 0.8413448

pnorm(1, mean=0, sd=1)

[1] 0.8413447

Differentiation

Examples

Density function = derivative of cumulative distribution function. Check with normal distribution.

library(numDeriv)
x <- seq(-2, 2, length.out=5)
grad(pnorm, x, mean=0, sd=1)

[1] 0.05399097 0.24197072 0.39894228 0.24197072 0.05399097

dnorm(x, mean=0, sd=1)

[1] 0.05399097 0.24197072 0.39894228 0.24197072 0.05399097

Derivative of exponential function is exponential function.

library(numDeriv)
exp(x)

[1] 0.1353353 0.3678794 1.0000000 2.7182818 7.3890561

grad(exp, x)

[1] 0.1353353 0.3678794 1.0000000 2.7182818 7.3890561

Delta method for standard error of function of parameter estimates

Maximum likelihood parameter estimate in Poisson regression

N  <- 100
X  <- rnorm(N, 0, 2)
mu <- exp(1 + 0.5*X)
Y  <- rpois(N, mu)
glmFit <- glm(Y ~ X, family=poisson(link="log"))
(bML   <- cbind(coef(glmFit)))

                 [,1]
(Intercept) 1.0960995
X           0.4842425

Function of the first two parameter estimates -> Evaluate function at ML estimate.

f <- function(b) { b[1] / b[2] }
f(bML)

[1] 2.263534

Delta method -> approximate standard error of function of parameter estimates -> approximate Wald confidence interval

library(numDeriv)
(bGrad <- grad(f, bML))

[1]  2.065081 -4.674382

(SEdelta <- sqrt(t(bGrad) %*% vcov(glmFit) %*% bGrad))

          [,1]
[1,] 0.2109704

(CIdelta <- c(lo=f(bML) - c(qnorm(1-(0.05/2))*SEdelta),
              up=f(bML) + c(qnorm(1-(0.05/2))*SEdelta)))

      lo       up 
1.850040 2.677029 

Standard errors from Hessian matrix

Poisson negative log-likelihood

nllPois <- function(b, X, Y) {
    mu <- exp(b[1] + b[2]*X)
    -sum(Y*log(mu) - mu)
}

Inverse of Hessian at ML estimate = variance-covariance matrix of parameter estimates.

h <- hessian(nllPois, x=bML, X=X, Y=Y)
solve(h)

              [,1]          [,2]
[1,]  0.0040165026 -0.0009308869
[2,] -0.0009308869  0.0004305852

Compare to vcov()

vcov(glmFit)

              (Intercept)             X
(Intercept)  0.0040165025 -0.0009308869
X           -0.0009308869  0.0004305852

Further resources

• CRAN Task View ‘Numerical Mathematics’

Useful packages

• statmod - gauss.quad() for Gaussian quadrature
• cubature
• polyCub

Detach (automatically) loaded packages (if possible)

try(detach(package:numDeriv))

Get the article source from GitHub

R markdown - markdown - R code - all posts