在R中通过卷积添加两个随机变量

Question

我想计算R中两个概率分布的卷积，并需要一些帮助。为了简单起见，假设我有一个变量x，其正态分布的平均值为1.0，标准差为0.5，以及一个log-normal分布的平均值为1.5，标准差为0.75的变量y。我想确定z = x + y。我知道z的分布事先是未知的。

顺便说一下，我正在处理的真实世界示例需要将两个随机变量相加，这两个随机变量按照许多不同的分布进行分配。

有人知道如何通过卷积x和y的概率密度函数来相加两个随机变量吗？

我尝试生成n个正态分布的随机值（使用上述参数）并将它们添加到n个log-normally分布的随机值中。但是，我希望知道是否可以使用卷积方法代替。任何帮助都将不胜感激。

编辑

谢谢你们的答案。我定义了一个PDF，并尝试进行卷积积分，但R在积分步骤上出现问题。我的PDF是Log Pearson 3，如下所示

dlp3 <- function(x, a, b, g) { 
 p1 <- 1/(x*abs(b) * gamma(a)) 
 p2 <- ((log(x)-g)/b)^(a-1) 
 p3 <- exp(-1* (log(x)-g) / b) 
 d <- p1 * p2 * p3 
 return(d) 
} 

f.m <- function(x) dlp3(x,3.2594,-0.18218,0.53441) 
f.s <- function(x) dlp3(x,9.5645,-0.07676,1.184) 

f.t <- function(z) integrate(function(x,z) f.s(z-x)*f.m(x),-Inf,Inf,z)$value 
f.t <- Vectorize(f.t) 
integrate(f.t, lower = 0, upper = 3.6)

R在最后一步抱怨，因为f.t函数是有界的，而我的积分限可能不正确。有没有什么解决办法？

- user3609848

---

我建议你查看distr包，或者至少快速浏览一下vignette。我认为这正是你要找的东西。虽然你用来生成随机值的策略也完全有效。 - MrFlick

非常感谢，@MrFlick。distr包是这样的任务的正确解决方案！ - API

Answer 1

这是一种方法。

f.X <- function(x) dnorm(x,1,0.5)        # normal (mu=1.5, sigma=0.5)
f.Y <- function(y) dlnorm(y,1.5, 0.75)   # log-normal (mu=1.5, sigma=0.75)
# convolution integral
f.Z <- function(z) integrate(function(x,z) f.Y(z-x)*f.X(x),-Inf,Inf,z)$value
f.Z <- Vectorize(f.Z)                    # need to vectorize the resulting fn.

set.seed(1)                              # for reproducible example
X <- rnorm(1000,1,0.5)
Y <- rlnorm(1000,1.5,0.75)
Z <- X + Y
# compare the methods
hist(Z,freq=F,breaks=50, xlim=c(0,30))
z <- seq(0,50,0.01)
lines(z,f.Z(z),lty=2,col="red")

使用distr包来完成相同的事情。

library(distr)
N <- Norm(mean=1, sd=0.5)          # N is signature for normal dist
L <- Lnorm(meanlog=1.5,sdlog=0.75) # same for log-normal
conv <- convpow(L+N,1)             # object of class AbscontDistribution
f.Z  <- d(conv)                    # distribution function

hist(Z,freq=F,breaks=50, xlim=c(0,30))
z <- seq(0,50,0.01)
lines(z,f.Z(z),lty=2,col="red")

- jlhoward

---

关于您的dlp3分布，您确定您的方程式正确吗？您的函数dlp3（...）发散。尝试x <- seq（0,100，.1）; plot（x，dlp3（x，3，-0.2，0.5））。实际上，可以证明它遵循大x时的~ x ^ 4 * log（x）^ 2。 - jlhoward

不确定您是否在5年后查看此内容，但（1）感谢您的点赞，这真的很有帮助。（2）您是否知道是否有一种方法可以反向操作？如果我有Z（或第二个示例中的f.Z）并且知道N（X），那么我能否将它们解卷积以获取L（Y）？我尝试了decon包，它有点作用，但是非常嘈杂。我似乎找不到任何分布内的东西。我想我可以手动完成，但似乎效果不佳 - 理想情况下 - 我需要使用具有变化sigma的N（X）进行操作 - 即异方差误差）。提前致谢。 - Mooks

Answer 2

我在使用不同密度参数时遇到了integrate()无法正常工作的问题，因此我想出了一种替代方法，使用黎曼逼近：

set.seed(1)

#densities to be convolved. could also put these in the function below
d1   <- function(x) dnorm(x,1,0.5) #
d2   <- function(y) dlnorm(y,1.5, 0.75) 

#Riemann approximation of convolution 
conv <- function(t, a, b, d) { #a to b needs to cover the range of densities above. d needs to be small for accurate approx.
  z <- NA
  x <- seq(a, b, d)
  for (i in 1:length(t)){
    print(i)
    z[i] <- sum(d1(x)*d2(t[i]-x)*d)
  }
  return(z)
}

#check against sampled convolution
X <- rnorm(1000, 1, 0.5)
Y <- rlnorm(1000, 1.5, 0.75)
Z <- X + Y

t <- seq(0, 50, 0.05) #range to evaluate t, smaller increment -> smoother curve

hist(Z, breaks = 50, freq = F, xlim = c(0,30))
lines(t, conv(t, -100, 100, 0.1), type = "s", col = "red")

- holdmygruyere