#' NNS Numerical Differentiation
#'
#' Determines numerical derivative of a given univariate function using projected secant lines on the y-axis. These projected points infer finite steps \code{h}, in the finite step method.
#'
#' @param f an expression or call or a formula with no lhs.
#' @param point numeric; Point to be evaluated for derivative of a given function \code{f}.
#' @param h numeric [0, ...]; Initial step for secant projection. Defaults to \code{(h = 0.1)}.
#' @param tol numeric; Sets the tolerance for the stopping condition of the inferred \code{h}. Defaults to \code{(tol = 1e-10)}.
#' @param digits numeric; Sets the number of digits specification of the output. Defaults to \code{(digits = 12)}.
#' @param print.trace logical; \code{FALSE} (default) Displays each iteration, lower y-intercept, upper y-intercept and inferred \code{h}.
#' @return Returns a matrix of values, intercepts, derivatives, inferred step sizes for multiple methods of estimation.
#' @author Fred Viole, OVVO Financial Systems
#' @references Viole, F. and Nawrocki, D. (2013) "Nonlinear Nonparametric Statistics: Using Partial Moments"
#' \url{https://www.amazon.com/dp/1490523995/ref=cm_sw_su_dp}
#' @examples
#' \dontrun{
#' f <- function(x) sin(x) / x
#' NNS.diff(f, 4.1)
#' }
#' @export
NNS.diff <- function(f, point, h = 0.1, tol = 1e-10, digits = 12, print.trace = FALSE){
Finite.step <- function(f, point, h){
f.x <- f(point)
f.x.h.min <- f(point - h)
f.x.h.pos <- f(point + h)
neg.step <- (f.x - f.x.h.min) / h
pos.step <- (f.x.h.pos - f.x) / h
return((c("f(x-h)" = neg.step,
"f(x+h)" = pos.step,
mean(c(neg.step, pos.step)))))
}
Bs <- numeric()
Bl <- numeric()
Bu <- numeric()
### Step 1 initalize the boundaries for B
### Initial step size h
f.x <- f(point)
f.x.h <- f(point - h)
### Y = mX + B
Y <- f.x
m <- (f.x - f.x.h)/h
mX <- ((f.x - f.x.h) / h) * point
B <- f.x - ((f.x - f.x.h) / h) * point
### Initial interval for B given inputted h-step-value
f.x.h.lower <- f(point - h)
f.x.h.upper <- f(point + h)
B1 <- f.x - ((f.x - f.x.h.lower) / h) * point
B2 <- f.x - ((f.x.h.upper-f.x) / h) * point
low.B <- min(c(B1, B2))
high.B <- max(c(B1, B2))
lower.B <- low.B
upper.B <- high.B
## Return "Derivative Does Not Exist" if lower.B and upper.B are identical to 20 digits
if(lower.B == upper.B){
original.par=par(no.readonly = TRUE)
par(mfrow = c(1, 2))
plot(f, xlim = c(point - (100 * h), point + (100 * h)), col = 'blue', ylab = 'f(x)')
points(point, f.x, pch = 19, col = 'red')
plot(f, xlim = c(point - 1, point + 1), col = 'blue', ylab = 'f(x)')
points(point, f.x, pch = 19, col = 'red')
par(original.par)
return(c("Derivative Does Not Exist"))}
new.B <- mean(c(lower.B, upper.B))
i <- 1L
while(i >= 1L){
Bl[i] <- lower.B
Bu[i] <- upper.B
new.f <- function(x) - f.x + ((f.x - f(point - x)) / x) * point + new.B
### SOLVE FOR h, we just need the negative or positive sign from the tested B
inferred.h <- uniroot(new.f, c(-2 * h, 2 * h), extendInt = 'yes')$root
if(print.trace) {print(c("Iteration" = as.integer(i), "h" = inferred.h, "Lower B" = lower.B, "Upper B" = upper.B))}
Bs[i] <- new.B
## Stop when the inferred h is within the tolerance level
if(abs(inferred.h) < tol) {
final.B <- mean(c(upper.B, lower.B))
slope <- solve(point, f.x - final.B)
z <- complex(real = point, imaginary = inferred.h)
original.par <- par(no.readonly = TRUE)
par(mfrow=c(1, 3))
## Plot #1
plot(f, xlim = c(min(c(point - (100 * h), point + (100 * h)), 0), max(c(point - (100 * h), point + (100 * h)), 0)), col = 'azure4', ylab = 'f(x)', lwd = 2, ylim = c(min(c(min(c(B1, B2)), min(na.omit(f((point - (100 * h)) : (point + (100 * h))))))), max(c(max(na.omit(f((point - (100 * h)) : (point + (100 * h))))), max(c(B1, B2))))), main = 'f(x) and initial y-intercept range')
abline(h = 0, v = 0, col = 'grey')
points(point, f.x, pch = 19, col = 'green')
points(point - h, f.x.h.lower, col = ifelse(B1 == high.B, 'blue', 'red'), pch = 19)
points(point + h, f.x.h.upper, col = ifelse(B1 == high.B, 'red', 'blue'), pch = 19)
points(x = rep(0, 2), y = c(B1, B2), col = c(ifelse(B1 == high.B, 'blue', 'red'), ifelse(B1 == high.B, 'red', 'blue')), pch = 1)
segments(0, B1, point - h, f.x.h.lower, col = ifelse(B1 == high.B, 'blue','red'), lty = 2)
segments(0, B2, point + h, f.x.h.upper, col = ifelse(B1 == high.B, 'red','blue'), lty = 2)
## Plot #2
plot(f, col = 'azure4', ylab = 'f(x)', lwd = 3, main = 'f(x) narrowed range and secant lines', xlim = c(min(c(point - h, point + h, 0)), max(c(point + h,point - h, 0))), ylim= c(min(c(B1, B2, f.x.h.lower, f.x.h.upper)), max(c(B1, B2, f.x.h.lower, f.x.h.upper))))
abline(h = 0, v = 0, col = 'grey')
points(point,f.x, pch = 19, col = 'red')
points(point - h, f.x.h.lower, col = ifelse(B1 == high.B, 'blue', 'red'), pch = 19)
points(point + h, f.x.h.upper, col = ifelse(B1 == high.B, 'red', 'blue'), pch = 19)
points(point, f.x, pch = 19, col = 'green')
segments(0, B1, point - h, f.x.h.lower, col = ifelse(B1 == high.B, 'blue', 'red'), lty = 2)
segments(0, B2, point + h, f.x.h.upper, col = ifelse(B1 == high.B, 'red', 'blue'), lty = 2)
points(x = rep(0, 2), y = c(B1, B2), col = c(ifelse(B1 == high.B, 'blue', 'red'), ifelse(B1 == high.B, 'red', 'blue')), pch = 1)
## Plot #3
plot(Bs, ylim = c(min(c(Bl, Bu)), max(c(Bl, Bu))), xlab = "Iterations", ylab = "y-inetercept", col = 'green', pch = 19, main = 'Iterated range of y-intercept')
points(Bl, col = 'red', ylab = '')
points(Bu, col = 'blue', ylab = '')
legend('topright', c("Upper y-intercept", "Lower y-intercept", "Mean y-intercept"), col = c('blue', 'red', 'green'), pch = c(1, 1, 19), bty = 'n')
par(original.par)
return(round(as.matrix(c("Value of f(x) at point" = f(point),
"Final y-intercept (B)" = final.B,
"DERIVATIVE" = slope,
"Inferred h" = inferred.h,
"iterations" = i,
Finite.step(f, point, h)[1 : 2],
"Averaged Finite Step Initial h " = Finite.step(f, point, h)[3],
"Inferred h" = Finite.step(f, point, inferred.h)[1 : 2],
"Inferred h Averaged Finite Step" = Finite.step(f, point, inferred.h)[3],
"Complex Step Derivative (Initial h)" = Im(f(z)) / Im(z))), digits))
}
## NARROW THE RANGE OF B BASED ON SIGN OF INFERRED.H
if(B1 == high.B){
if(sign(inferred.h) < 0) {
lower.B <- new.B
upper.B <- upper.B
} else {
upper.B <- new.B
lower.B <- lower.B
}
} else {
if(sign(inferred.h) < 0) {
lower.B <- lower.B
upper.B <- new.B
} else {
upper.B <- upper.B
lower.B <- new.B
}
}
new.B <- mean(c(lower.B, upper.B))
i <- i + 1
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.