# Negative value in density deconvolution

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))
}


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.

• It is about the decomposition of Round Trip Time(RTT) which was passively collected. The RTT is applied in TCP protocol, and could be affected by many factors. I want to decompose those effect, and build up a model with parameter to describe it. I inferred the main factor through samples and prior knowledge. Assumption that $X \sim Exp(\lambda)$ is a part of conditions, actually the model is far more complicated than that. I will provide that as complementary in the post above. @G.JayKerns – Readon Shaw Oct 5 '12 at 2:56
• @Readon Thanks for the background; the problem sounds interesting. I tried downloading your data (twice) from the link you gave julien, but it seems to be in some sort of binary format which I'm not sure how to read. If you'd export your data into a text file (or similar) I'd take a look at it. I'm not sure what you mean by "gap between the two basic distributions in the pdf", but if it means what it sounds like it means, then there might be something else you could do. – user1108 Oct 5 '12 at 12:10
• I am very glad that you and Julien are interested in this problem. I have convert the file from R workspace file to a zipped csv file, you can download it from here. The "gap" can be found in samples' histogram which can be used to estimate its pdf. @G.JayKerns – Readon Shaw Oct 5 '12 at 14:39

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

• 1.The estimator was not construct manually, I use density function provide as default in R. 2.the problem is related to density decomposition, there is no assumption about the "noise". 3. condition $X \sim Exp(\lambda)$ was mentioned above. And $Y$ as I guess (through the result with negative value), it might be mixture of two normal distribution.@julienstirnemann – Readon Shaw Oct 5 '12 at 3:00
• Sorry, still don't get it. From what I understand, your problem may be solved by nonparametric deconvolution estimators (for which R packages exist). Howevever, I'm not sure the conditions you suggest allow you to use those packages. Could you provide a minimal working example? Say a simple simulation program of what your data is and any prior conditions? – julien stirnemann Oct 5 '12 at 4:56
• I upload samples stored in R to a server here. The "value" elements in the variable blk3 is what I want to decomposed. There is no simulation program exist. @julienstirnemann – Readon Shaw Oct 5 '12 at 5:50
• Since this appears as a data analysis problem, you must be able to present your problem in a program format using simulation, in R for example, since you have been using that software. you may add your FFT code based on that example. – julien stirnemann Oct 5 '12 at 10:13
• I add up my code in above post. I don't know whether it is what you mean.@julienstirnemann – Readon Shaw Oct 5 '12 at 14:55