Is output of Deamer deconvolution not a density?

Question

I have a Model Y= X+e and need the density of X. The deamer package deconvolves the density for X, but if I use the simpsons rule to integrate this density, I get values which are above 1. The following example gives me a density which integral is 1.173454:

library(deamer) # deconvolution
library(Bolstad) # simpson rule

# The Y's I have are inv-Weibull distributed and the error's are inv-normal distributed.
# As the deconvolution of those would take a long time, i used runif in my example
# to simplify the problem. The following uncommented lines are what I like to deconvolve:

#library(actuar) # for rinvweibull
#y <- rinvweibull(30000, shape=5.53861156, scale=488)/1000
#y <- y[y<1.5]
#e <- 1/rnorm(30000, mean=0.0023853421, sd=0.0004784688)/1000
#e <- e[e<1.5]
#decon <- deamerSE(y, error=e, from=-0.1, to=0.3)

y <- runif(1000, min = 0.8, max=1.2)
e <- runif(1000, min = 0.1, max=0.5)
decon <- deamerSE(y, error=e, from=0.4, to=1)
plot(decon)

# following line gives me integral of density (with simpsons rule)
sintegral(decon$supp, decon$f)$value

I am not sure if this is just an estimation error or if I should consider downscaling the density so that the integrated density is 1:

# Downscaling
yValsScaled <- decon$f/area
    plot(decon$supp, yValsScaled, type="l")
(areaScaled <- sintegral(decon$supp, yValsScaled)$value)

What do you think?

btw: If you use a larger intervall with the from and to argments in deamerSE function, the integrated density will be even greater (because of the periodicity of the density). Usually I thought that with deamerSE I would get a density which integral (from -inf to inf) is approximately 1. Therefore I thougt that integrating the density using a smaller intervall (e.g. with from=0.4 and to=1 in the deamerSE function) should give me a density which integral is less than 1. But as you can see it doesn't. So I am rather confused.

Answer 1

It seems you are really putting my library at a test!! The output of deamer is an estimation of a density Several things:

1. Although, theoretically, any np density estimation (with deconvolution or not) ensures that the resulting estimate integrates to 1 (by use of a kernel which is itself a density), you have (necessary) numerical approximations at several levels in any algorithm (not to mention the numerical integration itself). Therefore, you cannot expect to find an integral exactly =1. Here 1.17 is reasonable given the sample size and the very weird problem you suggest.

2. You may test the deamer examples - just for checking - using integrate(predict, lower, upper, obj=est) for any deamer object named "est" and defining the "lower" and "upper" integration bounds. You will find that estimates integrate between 0.96 and 1.1 in my experience. By the way integrate() is built-in and quite reliable!

3. The problem you suggest is very awkward: deconvolving two uniforms with such measurement noise is very unusual and should probably be considered under parametric methods.

4. Your code does not add noise to your unobserved variable. Here is a corrected example, with a "reasonable" signal-to-noise ratio:

   x <- runif(2000, min = 0, max=1) 
   e <- runif(2000, min = - 0.1, max=0.1) 
   y <- x + e 
   err <- runif(1000, min = -0.1, max=0.1) 
   decon <- deamerSE(y, error=err, from=-1, to=2) 
   plot(decon)
   integrate(predict,-1,2,obj=decon)
   

This simple example actually yields an integral of 0.999...not too bad! Good luck with the rest! Julien Stirnemann