Search code examples
rdistributionfftnoiseconvolution

Deconvolution with R (decon and deamer package)


I have a model of the form: y = x + noise. I know the distribution of 'y' and of the noise and would like to have the distribution of 'x'. So I tried to deconvolve the distributions with R. I found 2 packages (decon and deamer) and I thought both methods should make more or less the same but I don't understand why deconvoluting with DeconPdf gives me a something like a normal distribution and deconvoluting with deamerKE gives me a uniform distribution. Here is an example code:

library(fitdistrplus) # for rweibull
library(decon) # for DeconPdf
library(deamer) # for deamerKE

set.seed(12345)
y <- rweibull(10000, shape=5.780094, scale=0.00204918)
noise <- rnorm(10000, mean=0.002385342, sd=0.0004784688)
sdnoise <- sd(noise)

est <- deamerKE(y, noise.type="Gaussian", 
                mean(noise), sigma=sdnoise)
plot(est)

estDecon <- DeconPdf(y, sdnoise, error="normal", fft=TRUE)
plot(estDecon)

Edit (in response to Julien Stirnemann):

I am not sure about re-parametrizing. My actual problem is: I have reaction time (RT) which theoretically can be described as f(RT) = g(discrimination time) + h(selection time), where f,g and h are can be transformations of those time values. I have "RT" and "discrimination time" values in my dataset. And I am interested in selection time or maybe h(selection time). With kernel density estimation I found out that the weibull distribution fits the 1/RT values best, while normal distribution fits 1/(discrimination time) best. That is why I can write my problem as 1/RT = 1/(discrimination time) + h(selection time) or y = x + noise (where I considered the noise to be 1/(discrimination time)). Simulating those reaction times gave me the following distribution with the following parameters:

y <- rweibull(10000, shape=5.780094, scale=0.00204918)
noise <- rnorm(10000, mean=0.002385342, sd=0.0004784688)

What do you mean with re-parametrizing? Using different values e.g. for the scale parameter?


Solution

  • As an answer to your last comment: the error shifts the observed values. The signal you wish to deconvolve is somewhere between 0 and ~ 0.3 i guess. Here is some code using deamer:

    library(actuar) # for rinvweibull
    library(deamer)
    set.seed(123)
    RT <- rinvweibull(30000, shape=5.53861156, scale=488)/1000
    RT <- RT[RT<1.5]
    noise <- 1/rnorm(30000, mean=0.0023853421, sd=0.0004784688)/1000
    noise <- noise[noise<1.5]
    
    ST <- deamerSE(RT, errors=noise, from=0, to=0.3)
    plot(ST)
    

    This is what you will get using non-parametric deconvolution (regardless of implementations, packages etc.). Just for your information, your signal-to-noise ratio is extremely low...what you observe is actually almost only noise. This strongly impacts on the estimation of the density you are interested in, especially when using nonparametric methods just like trying to find a needle in a haystack. You should reconsider estimating the density, and rather trying to get only a few quantities of interest...

    Good luck, Julien Stirnemann