generalized binary prior likelihood issue

[xxie6]

Hi! When using the generalized binary prior with flashier, I came across some instances where the KL estimate is negative. For example, for the flash fit generated by the code below:

set.seed(1)
X <- matrix(rnorm(100*1000), nrow=100, ncol=1000)
gb_flash_fit <- flash_init(X) |> flash_greedy(Kmax = 2, ebnm_fn = ebnm::ebnm_generalized_binary)

I found that all of the (negative) KL values were positive, e.g. one was 34.39.

I think there is an issue with the likelihood calculation. I looked through the code and I think in the generalized_binary.R file, line 228, the log(pi) term should be log(2*pi)?

[pcarbo]

@xxie6 The expression for the marginal likelihood is given in eq. 18 of the supplement of Yusha's paper. I must admit I don't understand how eq. 18 leads to line 228 — it looks like there are some intermediate calculations involved that I don't fully understand. You might want to try to write a function that takes as input "x", "s" and "fitted_g" (the output from calling the ebnm function), and apply eq. 18 to each pair (x[i], s[i]), and see if your calculation is the same as Yusha's. (While you are at it, you could compute other quantities such as the posterior mean as a way to reassure yourself.) I think it would be worthwhile — even if it ends up confirming that the calculation is correct, perhaps it would lead to some improvements in the generalized_binary.R code.