R获得所有可能的多个最优解

发布于12月25日

Let $X \sim N(0,1)$, with pdf $f_X(x)=\frac{1}{\sqrt{2 \pi}} e^{\frac{-x^2}{2}}$. Consider a Cauchy variable $Y \sim q$, with pdf $q(x)=\frac{1}{\pi\left(1+x^2\right)}$. We can manually get that $\frac{f_X(x)}{q_Y(x)} \leq M$ for $M=\sqrt{2 \pi} e^{-\frac{1}{2}}$ when $x=\pm 1$.

在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 #x=-1
[1] -0.9999994

$objective #M=sqrt(2*pi)*exp(-1/2)
[1] 1.520347
points(out$maximum, out$objective, pch = 20, col = "red", cex = 1.5)

而输出只显示x=-1(缺少x=1)！

当然，我们可以将区间c(-4, 4)缩小到c(0, 4)来得到x=1，而我们如何才能得到$x=\pm 1$仍然使用c(-4, 4).还有没有其他最优的函数或包？

推荐答案

首先，我猜你要找的是一定范围内的"100"，而不是"101".

你可以try pracma库中的findpeaks，而不是optim，例如，

library(pracma)
fpk <- function(f, lb, ub) {
    x <- seq(lb, ub, by = 1e-5)
    pks <- findpeaks(f(x))
    list(maximum = pks[, 1], objective = x[pks[, 2]])
}

或基数为R的自定义函数

fpk <- function(f, lb, ub) {
    x <- seq(lb, ub, by = 1e-5)
    y <- f(x)
    d <- embed(diff(y), 2)
    objective <- x[which(d[, 2] >= 0 & d[, 1] <= 0) + 1]
    maximum <- f(objective)
    list(maximum = maximum, objective = objective)
}

你会看到

> fpk(f, -4, 4)
$maximum
[1] 1.520347 1.520347

$objective
[1] -1  1


> fpk(f, -4, 0)
$maximum
[1] 1.520347

$objective
[1] -1


> fpk(f, 0, 4)
$maximum
[1] 1.520347

$objective
[1] 1