Masked density output for normal mixture

[nikhilmilind]

Bug Report

When I try to use the get_density function in the package to plot the density of a mixture of half-normal distributions, I get a (perfectly) symmetric distribution even though the input data is not symmetric around zero. The get_density function also cannot handle more than one input value when using the half-normal prior.

Replication Steps

First, I simulate some data that is clearly skewed.

set.seed(42)
b <- c(-rnorm(100, sd=5)^2, rnorm(100, sd=1)^2)
s <- runif(200, min=0.5)

Next, I fit ASH using the half-normal prior.

ash.res <- ash(b, s, mixcompdist="halfnormal")

I expect the following command to work, but it throws an error.

get_density(ash.res, seq(-10, 10, length.out=101))
# Error!

The function works for individual values, but the density is the same at -3 and 3 (very unlikely based on the data).

get_density(ash.res, -3)
get_density(ash.res, 3)

If I evaluate get_density one-value-at-a-time using sapply, I get a symmetric distribution.

x <- seq(-10, 10, length.out=101)
y <- sapply(x, function(i) { get_density(ash.res, i)$y })
plot(x, y, type="l")

Proposed Fix

I believe the issue lies in the comp_dens.tnormalmix function in R/tnormalmix.R file. Previously, this function returned a matrix of density values for each component (the matrix is component by target values). However, the output is not masked by the support of the mixture components. I add this masking in the fix.

Concerns

I'm not sure if this is the right place to implement the masking.

Replicating the Fix

First, generate the skewed data and fit the model.

set.seed(42)
b <- c(-rnorm(100, sd=5)^2, rnorm(100, sd=1)^2)
s <- runif(200, min=0.5)

ash.res <- ash(b, s, mixcompdist="halfnormal")

Now, get_density works properly with vector-valued targets.

get_density(ash.res, seq(-10, 10, length.out=101))

Furthermore, the density looks non-symmetric around zero.

density <- get_density(ash.res, seq(-10, 10, length.out=101))
plot(density$x, density$y, type="l")

[pcarbo]

Thank you @nikhilmilind for this feedback and for the fix. I will take a look.

I would also like to note that we have also developed the ebnm package which expands on many of the features of ashr, and has improved implementations of some of the different mixture priors. Perhaps you will end up finding ebnm more useful.

[pcarbo]

@nikhilmilind Your solution doesn't look quite right. But you are correct that there is a bug in this code.

I'll note that the "unimodal" distribution in the ebnm package should be more fully (and correctly) implemented, and might achieve a result similar to what you are looking for:

library(ebnm)
set.seed(42)
b <- c(-rnorm(100, sd=5)^2, rnorm(100, sd=1)^2)
s <- runif(200, min=0.5)
fit <- ebnm(b,s,prior_family = "unimodal")
plot(fit,incl_pm = FALSE,incl_cdf = TRUE)

I'll keep this PR open as I reminder to myself to fix this bug.

[nikhilmilind]

Excellent, thank you -- I'll try using ebnm in the meantime!