# Prediction Interval for Regression Models with Weights

I would like to understand how to manually derive a point-wise prediction confidence interval applicable to neural models.

Following http://cdn.intechopen.com/pdfs-wm/14915.pdf Confidence Intervals for Neural Networks and Applications to Modeling Engineering Materials,

my first aim was to construct a CI based on the Jacobian matrix related to model weights, using a simple lm model in R. By doing so I hope to get a comparable figure and then generalize it to non-linear (multilayer) set.

Below is my try. I am feeling really dumb in the very core of these calculations, and what I got did not fit with what I expected.

Question 1: is the Jacobian here relates to partial derivatives of the neural network output a rather then to a loss function of a?

Question 2: similarly, when it reads

the gradient of function a

does it relate to neural network output rather then it's loss?

rm(list=ls())

# load libraries ----------------------------------------------------------

list.of.packages <- c('data.table',
'numDeriv'
)
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(numDeriv)
library(data.table)

## create dummy data

nrows <- 1000

dat <- data.frame(x1 = rnorm(nrows, 0, 1)
, x2 = rnorm(nrows, 0, 1)
, x3 = rnorm(nrows, 0, 1)
)

dat$y = dat$x1 * 1 +
dat$x2 * 1 + dat$x3 * 1 +
rnorm(nrows, 0, 5)

lm_model <- lm(y ~ . -1, data = dat)

summary(lm_model)

## create function a (model output | W)

rm(fit_function)

fit_function <- function(x)
{

ww <- x

dat_copy <- as.data.table(copy(dat))

dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) ww[[x]])]

dat_copy[, y_hat := Reduce(+, Map(
function(x, y) x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]
)
)
]

# y <- mean(dat_copy[, (y - y_hat) ^ 2])

# X <- as.matrix(dat_copy[, y - y_hat])
#
# Xt <- t(X)
#
# XtX <- X %*% Xt
#
# sum_XtX <- sum(XtX) / 2

y <- unlist(dat_copy[, y_hat])

return(y)

}

## create function a (model output for one line | W)

rm(fit_function2)

fit_function2 <- function(x)
{

ww <- x

dat_copy <- as.data.table(copy(dat))

dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) ww[[x]])]

dat_copy[, y_hat := Reduce(+, Map(function(x, y)
x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]))]

y <- unlist(dat_copy[1, y_hat])

return(y)

}

Ww <- as.numeric(lm_model$coefficients) #I am using exactly fitted weights Jac_val <- numDeriv::jacobian( func = fit_function , x = Ww , method = "Richardson" , side = NULL ) T_Jac <- t(Jac_val) XtX <- T_Jac %*% Jac_val inv_X <- solve(XtX) #inversed matrix ## lets try derive a CI t_stats <- 3 dat_copy <- as.data.table(copy(dat)) dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) Ww[[x]])] dat_copy[, y_hat := Reduce(+, Map(function(x, y) x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]))] y_hat <- unlist(dat_copy[1, y_hat]) S <- sqrt(sum(dat_copy[, (y - y_hat) ^ 2]) / (nrow(dat) - (ncol(dat) - 1))) Grad_val <- numDeriv::grad( func = fit_function2 , x = Ww , method = "Richardson" , side = NULL ) StErr <- S * as.numeric(sqrt(1 + t(Grad_val) %*% inv_X %*% Grad_val)) preds <- stats::predict.lm( lm_model , newdata = dat[, 1:3] , se.fit = T , interval = "prediction" ) StErr_lm <- preds$se.fit[]

plot(preds$se.fit, type = 'l'); lines(preds$se.fit, type = 'p', col = 'red')


It's a clear mismatch in what I have done with what the lm SE is:

> StErr
 4.962455
> StErr_lm
 0.228429


## 1 Answer

I got it working, which is really cool. I was misunderstanding the se.fit part of the R::predict.lm. When I switched to interval = 'prediction' my math macthed the results.

To answer my own questions:

I need the partial derivatives of linear model function output with respect to given weight vector for both one observation modelled and all observations (to get the Jacobian matrix).

What I don't completely understand though is why

sign.level <- 0.05

stats::predict.lm(..., level = 1 - sign.level)

qt(p = 1 - sign.level / 2, df = nrow(dat) - 3, lower.tail = T, log.p = FALSE)


I don't have to divide 0.05 by 2 inside lm's function to get the same result as qt's in which I do have to divide 0.05 by 2 (for a two-tailed stats).

If matrix Jac_val %*% t(Jac_val) becomes singular and one cannot perform inversion, there is a workaround way to calculate solve(Jac_val %*% t(Jac_val)) that I also included in the code.

rm(list=ls())

# load libraries ----------------------------------------------------------

library('ggplot2')
library('numDeriv')
library('data.table')

## create dummy data

nrows <- 1000

dat <- data.frame(x1 = rnorm(nrows, 0, 1)
, x2 = rnorm(nrows, 0, 1)
, x3 = rnorm(nrows, 0, 1)
)

dat$y = dat$x1 * 1 +
dat$x2 * 1 + dat$x3 * 1 +
rnorm(nrows, 0, 5)

lm_model <- lm(y ~ . -1, data = dat)

summary(lm_model)

## create function a (model output | W)

rm(fit_function)

fit_function <- function(x)
{

ww <- x

dat_copy <- as.data.table(dat)

dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) ww[[x]])]

dat_copy[, y_hat := Reduce(+, Map(
function(x, y) x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]
)
)
]

y <- unlist(dat_copy[, y_hat])

return(y)

}

## create function a (model output for one line | W)

rm(fit_function2)

fit_function2 <- function(x)
{

ww <- x

dat_copy <- as.data.table(dat)

dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) ww[[x]])]

dat_copy[, y_hat := Reduce(+, Map(function(x, y)
x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]))]

y <- unlist(dat_copy[data_row, y_hat])

return(y)

}

Ww <- as.numeric(lm_model$coefficients) #I am using exactly fitted weights data_row <- 1 Jac_val <- numDeriv::jacobian( func = fit_function , x = Ww , method = "simple"# "Richardson" , side = NULL ) T_Jac <- t(Jac_val) XtX <- T_Jac %*% Jac_val inv_X <- solve(XtX) #inversed matrix ## alternative when XtX matrix is singular lambda <- 1e-5 nPara <- length(Ww) inv_X <- solve( t(Jac_val) %*% Jac_val + lambda * diag(nPara) ) %*% t(Jac_val) %*% Jac_val %*% solve( t(Jac_val) %*% Jac_val + lambda * diag(nPara) ) ## lets try derive a CI sign.level <- 0.05 dat_copy <- as.data.table(copy(dat)) dat_copy[, (paste0('w', seq(1, 3, 1))) := lapply(1:3, function(x) Ww[[x]])] dat_copy[, y_hat := Reduce(+, Map(function(x, y) x * y, .SD[, paste0('x', seq(1, 3, 1)), with = F], .SD[, paste0('w', seq(1, 3, 1)), with = F]))] y_hat <- unlist(dat_copy[data_row, y_hat]) S <- sqrt(sum(dat_copy[, (y - y_hat) ^ 2]) / (nrow(dat) - (ncol(dat) - 1))) Grad_val <- as.matrix( numDeriv::grad( func = fit_function2 , x = Ww , method = "simple"# "Richardson" , side = NULL ) ) StErr <- S * as.numeric(sqrt(1 + t(Grad_val) %*% inv_X %*% Grad_val)) y_hat_upp <- y_hat + StErr * qt(p = 1 - sign.level / 2, df = nrow(dat) - 3, lower.tail = T, log.p = FALSE) y_hat_low <- y_hat - StErr * qt(p = 1 - sign.level / 2, df = nrow(dat) - 3, lower.tail = T, log.p = FALSE) preds <- stats::predict.lm( lm_model , newdata = dat[, 1:3] #, se.fit = T , interval = "prediction" , level = 1 - sign.level ) y_hat_upp_lm <- preds[1, 3]  Visualizing: #dat_copy[, st_err := preds$se.fit]
dat_copy[, low_pred := preds[, 2]]
dat_copy[, upp_pred := preds[, 3]]
dat_copy[, step := 1:nrow(dat_copy)]

1 - sum(nrow(dat_copy[(y > upp_pred) | y < low_pred, ])) / nrow(dat_copy)

ggplot(data = dat_copy[step <= 500, ]) +

geom_line(aes(x = step, y = y_hat)
, size = 2
, alpha = 0.5
, color = 'red'
) +

geom_line(aes(x = step, y = y)
, size = 1
, alpha = 0.9
, color = 'black'
) +

geom_ribbon(aes(x = step, ymin = low_pred, ymax = upp_pred)
, color = 'red'
, fill = 'red'
, alpha = 0.1
, size = 0.2
) +

theme_bw()