识别曲线的全局最小值和最大值

Question

我有下面的函数

GetExtrema

（感谢@ZheyuanLi），我用它来识别曲线最小值和最大值处的x（年龄）和y（BMI，kg/m2）。

因为可能存在局部（而不是全局）最小值/最大值，我如何编辑

GetExtrema

函数以要求曲线的一阶导数之前0.25年（即3个月）为正，0.25年为负在达到最大值后，它被视为真正的最大值（最小值也相同）。

# Load packages

library(tidyverse)
library(ggeffects)
library(mgcv)

# data

dat <- read.csv((
  "https://raw.githubusercontent.com/aelhak/NCRM2023/main/bmi_long.csv")) %>% 
  select(age, bmi)

# model

model <- gam(
  bmi ~ s(sqrt(age), bs = 'ps'), 
  method = "REML", data = dat
  )

# predictions

pred <- predict_response(model, terms = "age [all]") %>% 
  as.data.frame() %>% select(x, predicted) %>% 
  rename(y = predicted)

with(pred, plot(x, y, type = "l"))

# Identify x and y at curve min and max

GetExtrema <- function (x, y) {
  turn <- diff.default(sign(diff.default(y)))
  maxInd <- which(turn == -2) + 1L
  minInd <- which(turn == 2) + 1L
  list(min.x = x[minInd], min.y = y[minInd],
       max.x = x[maxInd], max.y = y[maxInd])
}

GetExtrema(pred$x, pred$y)

Answer 1

首先，你不是在寻找全局极值。显然，全球最低年龄是最低年龄。您正在寻找具有额外条件的局部极值。

gratia 包提供了一个使用有限差分方法获取导数的函数。我将使用这个函数来实现极值搜索。然而，数字的本质是您需要提供容差，而这些容差将始终取决于数据（此处为模型预测）。

library(gratia)

GetExtrema2 <- function(model, x, tolx, tolmax, tolmin, eps = 1e-07) {
  x <- seq(min(x), max(x), by = tolx)
  
  #dy/dx
  d1 <- derivatives(model, 
                    data = data.frame(age = x), 
                    order = 1, type = "central", eps = eps)
  
  #d^2y/dx^2
  d2 <- derivatives(model, 
                    data = data.frame(age = x), 
                    order = 2, type = "central", eps = eps)
  
  d2chk <- round(d2$.derivative, -log10(tolmax)) < 0
  xmax <- (d1$`sqrt(age)`^2)[abs(d1$.derivative) < tolmax & d2chk]
  
  d2chk <- round(d2$.derivative, -log10(tolmin)) > 0
  xmin <- (d1$`sqrt(age)`^2)[abs(d1$.derivative) < tolmin & d2chk]
  
  ymax <- unname(predict(model, newdata = data.frame(age = xmax)))
  ymin <- unname(predict(model, newdata = data.frame(age = xmin)))
  
  plot(x, predict(model, newdata = data.frame(age = x)), type = "l")
  points(xmin, ymin, col = "blue", pch = 16)
  points(xmax, ymax, col = "dark red", pch = 16)
  
  list(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax)
}

print(ex <- GetExtrema2(model, x = dat$age, tolx = 1e-2, tolmax = 1e-2, tolmin = 1e-3))
#$xmin
#[1] 4.626306 4.636306
#
#$ymin
#[1] 15.69887 15.69887
#
#$xmax
#[1] 0.8263063
#
#$ymax
#[1] 17.21079

如您所见，根据提供的公差，该函数找到一个最大点和两个最小点（y 值几乎相同）。您可以调整容差以找到唯一的最小值，或者将结果四舍五入并使用

unique

。只需查看剩余极值之间的距离即可实现您的条件。但是，如果您得到“过于局部”的极值，您可能只想减小平滑基础的维度（请参阅

k

的

s

参数）。