Calculate the Derivative of a Data-Function in R

Calculate the derivative of a data-function in r

Was also going to suggest an example of a smoothed spline fit followed by prediction of the derivative. In this case, the results are very similar to the diff calculation described by @dbaupp:

spl <- smooth.spline(x, y=ycs)
pred <- predict(spl)

plot (x, ycs, log="xy")
lines(pred, col=2)

ycs.prime <- diff(ycs)/diff(x)
pred.prime <- predict(spl, deriv=1)

plot(ycs.prime)
lines(pred.prime$y, col=2)

R Function to Find Derivative of Every Point in Time Series

Your proposed example using the BJSales data is decidedly not differentiable,

so instead I will show the derivative of a much smoother function. If your real data is smooth, this should work for you.

The simplest way to approximate the derivative is simply to use finite differences.

f'(x) ≈ (f(x+h) - f(x))/h

## Smooth sample function
x = seq(0,10,0.1)
y = x/2 + sin(x)
plot(x,y, pch=20)

## Simplest - first difference
d1 = diff(y)/diff(x)
d1 = c(d1[1],d1)

Let's use it to plot a tangent line as an error check. I picked a place to draw the tangent line arbitrarily: the 18th point, x=1.7

plot(x,y, type="l")
abline(y[18]-x[18]*d1[18], d1[18]) 

To get the data.frame that you requested, you just need

Derivative = data.frame(x, d1)

get a derivative by knowing two numeric vectors

To find the derivative use the numeric approximation: (y2-y1)/(x2-x1) or dy/dx. In R use the diff function to calculate the difference between 2 consecutive points:

x<-rnorm(100)
y<-x^2+x

#find the average x between 2 points
avex<-x[-1]-diff(x)/2
#find the numerical approximation
#delta-y/delta-x
dydx<-diff(y)/diff(x)

#plot numeric approxiamtion
plot(x=avex, dydx)
#plot analytical answer
lines(x=avex, y=2*avex+1)

Calculate Derivative of Function at a Point in R

Just install the package numDeriv and use the grad function. Here are a few simple examples that are easy to check.

library(numDeriv)

grad(sin, 1:3)
[1]  0.5403023 -0.4161468 -0.9899925
cos(1:3)
[1]  0.5403023 -0.4161468 -0.9899925

f = function(x) x^2 + 2*x +3
grad(f, 1:3)
[1] 4 6 8
2*(1:3) + 2
[1] 4 6 8

R - how to calculate derivatives of two-dimensional discrete function?

So, finally I solved it (not the best solution, but it could help someone anyway). I decided to start with numerical differentiation (central differences) and not to use local spline/polynomial fitting as @42- suggested.

# our dataset consists of points (x, y, z), where z = f(x, y)

z_dx <- function(x, y) { # first partial derivative dz/dx
    pos_x <- which(df_x == x)

    if (pos_x == 1 | pos_x == length(df_x)) return(NA)

    xf <- df_x[pos_x + 1]
    xb <- df_x[pos_x - 1]

    return((df[which(df$x == xf & df$y == y), "z"] - df[which(df$x == xb & df$y == y), "z"]) / (xf - xb))
}

z_dy <- function(x, y) { # first partial derivative dz/dy
    pos_y <- which(df_y == y)

    if (pos_y == 1 | pos_y == length(df_y)) return(NA)

    yf <- df_y[pos_y + 1]
    yb <- df_y[pos_y - 1]

    return((df[which(df$y == yf & df$x == x), "z"] - df[which(df$y == yb & df$x == x), "z"]) / (yf - yb))
}

z_dx2 <- function(x, y) { # second partial derivative d2z/dx2
    pos_x <- which(df_x == x)

    if (pos_x <= 2 | pos_x >= (length(df_x) - 1)) return(NA)

    xf <- df_x[pos_x + 1]
    xb <- df_x[pos_x - 1]

    return((df[which(df$x == xf & df$y == y), "dx"] - df[which(df$x == xb & df$y == y), "dx"]) / (xf - xb))
}

z_dy2 <- function(x, y) { # second partial derivative d2z/dy2
    pos_y <- which(df_y == y)

    if (pos_y <= 2 | pos_y >= (length(df_y) - 1)) return(NA)

    yf <- df_y[pos_y + 1]
    yb <- df_y[pos_y - 1]

    return((df[which(df$y == yf & df$x == x), "dy"] - df[which(df$y == yb & df$x == x), "dy"]) / (yf - yb))
}

z_dxdy <- function(x, y) { # second partial derivative d2z/dxdy
    pos_y <- which(df_y == y)

    if (pos_y <= 2 | pos_y >= (length(df_y) - 1)) return(NA)

    yf <- df_y[pos_y + 1]
    yb <- df_y[pos_y - 1]

    return((df[which(df$y == yf & df$x == x), "dx"] - df[which(df$y == yb & df$x == x), "dx"]) / (yf - yb))
}

# read dataframe and set proper column names
df <- read.table("map.dat", header = T)
names(df) <- c("x", "y", "z")

# extract unique abscissae
df_x <- unique(df$x)
df_y <- unique(df$y)

# calculate first derivatives numerically
df$dx <- apply(df, 1, function(x) z_dx(x["x"], x["y"]))
df$dy <- apply(df, 1, function(x) z_dy(x["x"], x["y"]))

# set tolerance limits, so we can filter out every non-stationary point later
tolerance_x <- 30.0
tolerance_y <- 10.0

# select only stationary points
df$stat <- apply(df, 1, function(x) ifelse(is.na(x["dx"]) | is.na(x["dy"]), FALSE, ifelse(abs(x["dx"]) <= tolerance_x & abs(x["dy"]) <= tolerance_y, TRUE, FALSE)))
stationaries <- df[which(df$stat == TRUE), ]

# compute second derivatives for selected points
stationaries$dx2 <- apply(stationaries, 1, function(x) z_dx2(x["x"], x["y"]))
stationaries$dy2 <- apply(stationaries, 1, function(x) z_dy2(x["x"], x["y"]))
stationaries$dxdy <- apply(stationaries, 1, function(x) z_dxdy(x["x"], x["y"]))

# let's make classifying function...
stat_type <- function(dx2, dy2, dxdy) {
    if (is.na(dx2) | is.na(dy2) | is.na(dxdy)) return("unk")
    if (dx2 * dy2 - dxdy^2 < 0) return("saddle")
    if (dx2 * dy2 - dxdy^2 > 0) {
        if (dx2 < 0 & dy2 < 0) return("max")
        else if (dx2 > 0 & dy2 > 0) return("min")
        else return("unk")
    }
    else return("unk")
}

# ...and apply it to our stationary points
stationaries$type <- apply(stationaries, 1, function(x) stat_type(x["dx2"], x["dy2"], x["dxdy"]))
stationaries <- stationaries[which(stationaries$type != "unk"), ] # leave out all non-stationarities

# set upper cutoff (in the units of z = f(x, y)) - points with the greatest value of z can form very large areas where derivatives are so close to zero that we can't discriminate them (this is often the case when our map is built by numerical integration with predefined grid so final result is discrete function of two variables)
cutwidth <- 2.5
stationaries <- stationaries[which(stationaries$z <= (max(stationaries$z) - cutwidth)), ] # select only points which are lying below maximum by selected cutoff

# create dataframes for minima, maxima and saddles
minima <- stationaries[which(stationaries$type == "min"), ]
maxima <- stationaries[which(stationaries$type == "max"), ]
saddles <- stationaries[which(stationaries$type == "saddle"), ]

Things becomes easier if your map have constant grid size - then there is no need to calculate stepsize for differentiation for every point

---

Related Topics