R - 如何计算二维离散函数的导数？
R - how to calculate derivatives of two-dimensional discrete function?
我有二维离散函数——实际上是一个用两个变量定义的 3D 表面。这个表面是由另一个程序构建的，因此我的数据文件中有大量 (X, Y, Z) 点。它可以用 ggplot 成功绘制，这很好。
但现在我应该找到这些曲面上的所有固定点 - 最大值、最小值和鞍点。我知道这样的任务可以借助一阶和二阶偏导数准则轻松解决，但我不知道如何在实践中实现它。
至于我，我认为最明显的方法是在两个方向上穿过我的表面，计算数值导数，然后测试获得的值是否为零（具有一定的公差）。我找到了 R 的 numDeriv 包，但它似乎希望 func 对象作为实值函数，而我只有离散的点集。
所以我想知道是否有任何其他建议、更好的方法或者甚至可能是现成的包来解决这样的任务？谢谢！
在 Rhelp 上搜索，我发现了一个类似的问题（虽然不是很广泛）。它只是要求找到局部最大值，但推广到其他情况应该不会太困难。提问者提供了一个数据集：
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- 建议的那样使用局部 spline/polynomial 拟合。
# 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"), ]
如果您的地图具有恒定的网格大小，事情就会变得更容易 - 然后就不需要为每个点的差异计算步长
我有二维离散函数——实际上是一个用两个变量定义的 3D 表面。这个表面是由另一个程序构建的，因此我的数据文件中有大量 (X, Y, Z) 点。它可以用 ggplot 成功绘制，这很好。
但现在我应该找到这些曲面上的所有固定点 - 最大值、最小值和鞍点。我知道这样的任务可以借助一阶和二阶偏导数准则轻松解决，但我不知道如何在实践中实现它。
至于我，我认为最明显的方法是在两个方向上穿过我的表面，计算数值导数，然后测试获得的值是否为零（具有一定的公差）。我找到了 R 的 numDeriv 包，但它似乎希望 func 对象作为实值函数，而我只有离散的点集。
所以我想知道是否有任何其他建议、更好的方法或者甚至可能是现成的包来解决这样的任务？谢谢！
在 Rhelp 上搜索，我发现了一个类似的问题（虽然不是很广泛）。它只是要求找到局部最大值，但推广到其他情况应该不会太困难。提问者提供了一个数据集：
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- 建议的那样使用局部 spline/polynomial 拟合。
# 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"), ]
如果您的地图具有恒定的网格大小，事情就会变得更容易 - 然后就不需要为每个点的差异计算步长