用R解卷积(decon和deamer程序包) [英] Deconvolution with R (decon and deamer package)

问题描述

我有一个形式的模型:y = x +噪声.我知道'y'和噪声的分布，并且想知道'x'的分布.因此，我尝试用R对卷积进行反卷积.我发现了2个软件包(decon和deamer)，我认为这两种方法应大致相同，但我不明白为什么用DeconPdf反卷积会给我类似正态分布的分布，用deamerKE去卷积可以得到均匀的分布.这是示例代码:

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)

---

编辑(针对 Julien Stirnemann ):

---

Edit (in response to Julien Stirnemann):

我不确定是否要重新设置参数.我的实际问题是: 我有反应时间(RT)，理论上可以描述为f(RT)= g(判别时间)+ h(选择时间)，其中f，g和h可以是这些时间值的变换. 我的数据集中有"RT"和歧视时间"值.而且我对选择时间或h(选择时间)感兴趣.通过核密度估计，我发现威布尔分布最适合1/RT值，而正态分布最适合1/(区分时间). 这就是为什么我可以这样写我的问题的原因:1/RT = 1/(判别时间)+ h(选择时间)或y = x +噪声(这里我认为噪声为1/(判别时间)).通过模拟这些反应时间，可以得到具有以下参数的以下分布:

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?

推荐答案

作为对您的最后评论的答案:错误会改变观察值.我想您希望解卷积的信号在0到〜0.3之间.这是一些使用deamer的代码:

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...

祝你好运， 朱利安·斯特尼曼(Julien Stirnemann)

Good luck, Julien Stirnemann

这篇关于用R解卷积(decon和deamer程序包)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！