# 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

## Detach (automatically) loaded packages (if possible)

try(detach(package:numDeriv))