Likelihood not maximized using prior_family = "normal"

[jean997]

I have found some cases where the fitted_g returned by ebnm does not maximize the log likelihood when using prior_family = "normal". For example, in the code below, ebnm finds the sd to be 0.107 but it should be 0.193.

set.seed(1)
n <- 1e5
A <- 0.03
b_hat <- rnorm(n =n, mean = 0, sd = sqrt(1 + A))
b_post <- ebnm(x = b_hat, s = 1, prior_family = "normal")
b_post$fitted_g$sd
# [1] 0.1066113
b_post$log_likelihood
'log Lik.' -143739.4 (df=1)

sigma2_max <- mean(b_hat^2) - 1
sqrt(sigma2_max)
# [1] 0.1930426
sum(dnorm(b_hat, mean = 0, sd = sqrt(1 + sigma2_max), log = TRUE))
# [1] -143723.2

This doesn't seem to always be an issue. For example if A = 0.3 instead of 0.03, the answers agree out to 6 decimal points.

set.seed(1)
n <- 1e5
A <- 0.3
b_hat <- rnorm(n =n, mean = 0, sd = sqrt(1 + A))
b_post <- ebnm(x = b_hat, s = 1, prior_family = "normal")
b_post$fitted_g$sd
# [1] 0.5560356
b_post$log_likelihood
'log Lik.' -155363.5 (df=1)

sigma2_max <- mean(b_hat^2) - 1
sqrt(sigma2_max)
# [1] 0.5560306
sum(dnorm(b_hat, mean = 0, sd = sqrt(1 + sigma2_max), log = TRUE))
[1] -155363.5

[pcarbo]

@jean997 It appears that nlm without the Hessian fails — I'm not quite sure why. Perhaps Jason has an idea.

library(ebnm)
set.seed(1)
n <- 100000
A <- 0.03
b_hat <- rnorm(n = n,mean = 0,sd = sqrt(1 + A))
b_post <- ebnm(x = b_hat,s = 1,prior_family = "normal",
                optmethod = "nohess_nlm")
b_post$fitted_g$sd
# [1] 0.1066
b_post <- ebnm(x = b_hat,s = 1,prior_family = "normal",
                optmethod = "nlm")
b_post$fitted_g$sd
# [1] 0.193

[jean997]

Thanks!