R - как вычислить производные двумерной дискретной функции?

В поисках Релпа я нашел похожий вопрос (хотя и не такой широкий). Он просто попросил найти локальные максимумы, но обобщение на другие случаи не должно быть слишком сложным. Опрашивающий предоставил набор данных:

require(MASS)
topo.lo <- loess(z ~ x * y, topo, degree = 1, span = 0.25,
                 normalize = FALSE)
topo.mar <- list(x = seq(0, 6.5, 0.1), y = seq(0, 6.5, 0.1))
new.dat <- expand.grid(topo.mar)
topo.pred <- predict(topo.lo, new.dat)
## draw the contour map based on loess predictions


library(rgl)


persp3d(topo.mar$x, topo.mar$y, topo.pred, shade=0.5, col="blue")

И получил этот ответ (который, по общему признанию, не рассчитывает производные, а числовые производные являются лишь первыми отличиями), что и используется в этом тесте:

hasmax <- function(mtx, x, y)  if(  (mtx[x,y] > mtx[x,y-1]) &
                                    (mtx[x,y] > mtx[x,y+1]) &
                                    (mtx[x,y] > mtx[x-1,y]) &
                         (mtx[x,y] > mtx[x+1,y]) ) {return(TRUE ) } else {return(FALSE)}


for(i in 3:(dim(topo.pred)[1] -4)) {
    for(j in 3:(dim(topo.pred)[2]-4) )  {
        if( hasmax(topo.pred, i , j) ){print(c(topo.mar$x[i],topo.mar$y[j]))}  }}

Я публикую его с небольшим изменением, потому что это был мой ответ: -)

Я думаю, что могут быть пакеты, которые имеют дело с пространственными наборами данных, но я ' Я не особенно знаком, поэтому вы можете искать в представлении задач CRAN.

Итак, наконец-то я решил это (не лучшее решение, но оно все равно может кому-то помочь). Я решил начать с численного дифференцирования (центральных различий) и не использовать локальную сплайн / полиномиальную подгонку, как предлагалось в @ 42- .

# 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"), ]

Все становится проще, если ваша карта имеет постоянный размер сетки - тогда нет необходимости рассчитывать размер шага для дифференциации для каждой точки