Negative value in density deconvolution

Question

I have a set of samples from measurement. It can be expressed as random variable $Z = X + Y$. So that $Z$ is observable, but $X,Y$ is unobservable. According to prior knowledge and histogram of $Z$, I can assume that $X$ following exponential distribution with a manual selected parameter. The estimation of probability density function(pdf) of $Y$ is what I wanted. I use fft method to process this density deconvolution problem. More details can be found in this post.

To solve this numerical density deconvolution problem we adopt "deconv" function supplied by pracma package in R. The results contain a block of negative value which conflict the definition of probability density. If I ignore all of those invalid value by setting them to zero and convolute it back, there would be a great mount of differences between reconvoluted value and samples. Why? How to process it?

Furthermore, if I want to get parameters of $Y$ after its pdf was estimated, perhaps a mixture of normal distribution, what should I do?

=================================================

The background is provided as complementary: The samples are collected Round Trip Time(RTT) marked as $D_{rtt}$. It can be expressed as $D_{rtt} = D_{mech} + D_{delay}$, where $D_{mech}$ refers to time delay caused by mechanism, $D_{delay}$ is the real delay which is the main task for this analysis. $D_{mech}$ can be concluded as mixture distribution of exponential distribution and unifrom distribution. Fortunately, there is a gap between the two basic distributions in the pdf. So that I can process it separately as I modeled above.

=================================================

I add up my code for sample processing in R below.

delay_est = function(x, samples = 1000, thresh = 0.1){  
    library(pracma)
    lp = subset(x, value < thresh)
    d = hist(lp$value, samples)
len = length(d$density)

    fit = function(lambda){
      di = dexp(d$mids, lambda)
  dx = deconv(c(d$density, rep(0, len)), subset(di, di > 0))
      rd = dx$q[dx$q<0]
      return(sum(abs(rd)))
    }

    opt = optim(500, fit, method='Brent', lower=0, upper=10000)
    print(str(opt))  
    di = dexp(d$mids, opt$par)      
    dx = deconv(c(d$density, rep(0, len)), di)
return(list(dx$q, opt))
  }

Answer 1

The short answer to your question is, "Maybe it's not possible." I mean that, depending on $Z$, it might be impossible to decompose it into $X + Y$ such that $X$ and $Y$ are independent, $X$ is exponential, and $Y$ is (a valid) something else.

Part of my research years ago centered on a related but more specific question: given an ordinary $Z$, is it possible to deconvolve into $X + Y$ where $X$ and $Y$ are independent and identically distributed (IID)? It turns out, depending on $Z$, the answer is "No." Well, that's only half true: The answer is "Yes", but sometimes the distribution of $X$ (and $Y$) is negative.

I strongly suspect that a similar thing is happening to you. You're given $Z$, yet you are specifying $X$ to be exponential and independent of $Y$. As you have discovered, the numerical estimation of the PDF of $Y$ turns negative. That means there is something about $Z$ which is incompatible with your other two requirements while maintaining that $Y$ is a random variable in the classical sense of the term.

Deconvolution is a difficult problem. It's been a few years; I would need to go back and brush up on what I knew. Also, I'm not really current on developments in recent years (if any).

But to be clear: there is nothing mathematically inconsistent or illegal about an estimated function turning negative on some parts of its domain. There hasn't been a lot of work done concerning statistical inference in cases like that, so you're on your own when it comes time to estimate parameters, fit a model, etc.

I'd be interested to hear more about the physical problem/application that gave rise to this question: details, data, and so forth. Who knows, maybe you've come up with something new and really cool.

Answer 2

It is possible to get negative values depending on the projection basis of functions your are considering when using fft methods in density deconvolution problems. However: could you give more details regarding your problem? how have you constructed your estimator? is your problem related to denoising? Does the law of $X$ belong to a knwon parametric class?
cheers