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$.
In R, I have tried
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)
while the output only shows x=-1 (missing x=1)!
Surely, we can shrink the interval c(-4, 4) to c(0, 4) to get x=1, while how we can get both $x=\pm 1$ still using c(-4, 4). Any other optimal functions or packages?
First of all, I guess what you are after is "local maximums" within a given range, instead of "global maximum".
You can try findpeaks from pracma library, instead of optim, e.g.,
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]])
}
or a custom function with base 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)
}
and you will see
> 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