如何在R语言中获取f(x)=dnorm(x)/dcauchy(x)的全部最优解

问题背景 #

设X服从标准正态分布N(0,1)，概率密度函数为$f_X(x)=\frac{1}{\sqrt{2 \pi}} e^{\frac{-x2}{2}}$；Y服从柯西分布，概率密度函数为$q(x)=\frac{1}{\pi\left(1+x^2\right)}$。经手动计算可知，$\frac{f_X(x)}{q_Y(x)}$在$x=\pm1$时取得最大值$M=\sqrt{2 \pi} e^{-\frac{1}{2}}$。

执行以下R代码时：

f <- function(x) dnorm(x) / dcauchy(x)
curve(f, -4, 4, n = 200, col = 4); grid()
out = optimize(f, interval = c(-4, 4), maximum = TRUE)
out
# $maximum
# [1] -0.9999994
# 
# $objective
# [1] 1.520347
points(out$maximum, out$objective, pch = 20, col = "red", cex = 1.5)

仅得到$x=-1$这个最大值点，缺失$x=1$。虽然缩小区间到c(0,4)可得到$x=1$，但需要在区间c(-4,4)下同时获取$x=\pm1$两个最优解。

解决方案 #

方法1：拆分区间分别搜索 #

由于目标函数是偶函数（$f(-x)=f(x)$），可以将原区间拆分为左右两部分分别搜索最大值，再合并结果：

f <- function(x) dnorm(x) / dcauchy(x)

# 左半区间搜索
left_opt <- optimize(f, interval = c(-4, 0), maximum = TRUE)
# 右半区间搜索
right_opt <- optimize(f, interval = c(0, 4), maximum = TRUE)

# 查看两个最优解
list(left_max = left_opt, right_max = right_opt)

执行后会同时得到$x\approx-1$和$x\approx1$两个最大值点，以及对应的最大值。

方法2：使用多极值点检测工具包 #

可以借助pracma包的findpeaks函数先定位离散点中的峰值位置，再对每个峰值附近做精细优化，适合处理更复杂的多极值函数：

library(pracma)

f <- function(x) dnorm(x) / dcauchy(x)
# 生成足够密集的离散点
x_grid <- seq(-4, 4, length.out = 1000)
y_vals <- f(x_grid)

# 检测峰值位置（返回包含峰值x坐标的矩阵）
peaks <- findpeaks(y_vals, x = x_grid, sortstr = TRUE)

# 对每个峰值附近做精细优化
optimal_results <- lapply(peaks[, 1], function(peak_x) {
  optimize(f, interval = c(peak_x - 0.5, peak_x + 0.5), maximum = TRUE)
})

optimal_results

为什么optimize无法直接返回多个极值？ #

R内置的optimize()函数基于Brent单变量优化算法，它只能找到区间内的一个局部极值点（当存在多个极值时，收敛结果取决于算法的初始搜索方向），因此无法直接返回所有极值点。

内容的提问来源于stack exchange，提问作者John Stone