When mu is fixed, check the gradient at pi0 = 0 and pi0 = 1. This would have the nice side effect of fixing a minimum scale for the slab, which would help stabilize the problem for pi0 near 1.
pn_llik_fn = function(x) {
if (x[1] >= 1) {
return(sum(dnorm(dat, sd = s, log = TRUE)))
} else {
return(sum(
logscale_add(log(x[1]) + dnorm(dat, sd = s, log = TRUE),
log(1 - x[1]) + dnorm(dat, sd = sqrt(s^2 + x[2]), log = TRUE))
))
}
}
numDeriv::grad(pn_llik_fn, x = c(1, min_sigma2))