基 R 函数卷积

Question

我正在构建一个函数集合，这些函数从两个独立随机变量的 pdf 中返回概率密度函数 (pdf)。

最常见的例子是独立随机变量 X、Y 的总和，由它们的 pdf 卷积给出。

在post 之后，我定义了以下函数，该函数将一对 pdf 作为参数，并返回它们的卷积：

dSumXY <- function(dX, dY){

  # Create convolution of distributions.
  dReturn <- function(z){
    integrate( function(x,z){
      dX(x) * dY(z - x) },
      -Inf, Inf, z)$value
  }

  # Vectorize convolution.
  dReturn <- Vectorize(dReturn)

  return(dReturn)  
}

这在以下示例中按预期工作：

# Define pdfs of two (identical) uniform [-1,1] distributions
unifX <- function(x) dunif(x, min = -1, max = 1)
unifY <- function(x) dunif(x, min = -1, max = 1)

# Find convolution of their pdfs.
convXY <- dSumXY(unifX, unifY)

# Plot the convolved pdf.
plot(seq(-3,3,by = 0.1), convXY(seq(-3,3,by = 0.1) ), type = 'l')

# Sample from the distribution
convSample <- runif(10000, min = -1, max = 1) + runif(10000, min = -1, max = 1)

# Plot density of sample.
lines( density(convSample) , col = "red" )

更一般地说，这适用于许多均匀分布对的组合，但是当我尝试对一对 Uniform[1,2] 分布进行卷积时，我没有得到真正的结果：

# Define pdfs of two (identical) uniform [1,2] distributions
unifX2 <- function(x) dunif(x, min = 1, max = 2)
unifY2 <- function(x) dunif(x, min = 1, max = 2)

# Find convolution of their pdfs.
convXY2 <- dSumXY(unifX2, unifY2)

# Plot the convolved pdf.
plot(seq(1,5,by = 0.1), convXY2(seq(1,5,by = 0.1) ), type = 'l')

# Sample from the distribution
convSample2 <- runif(10000, min = 1, max = 2) + runif(10000, min = 1, max = 2)

# Plot density of sample.
lines( density(convSample2) , col = "red" )

特别是，（有一点概率知识！）很明显，两个 Uniform[1,2] 变量之和的 pdf 在 3.75 处应为非零；特别是它应该等于 1/4。不过

# This should be equal to 1/4, but returns 0.
convXY2(3.75)

我已经在两台不同的机器上尝试过这个并复制了同样的问题，所以我很想看看问题是从哪里引起的。

提前致谢。

Answer 1

问题来自 Integrate 函数，该函数正在努力解决该函数具有紧凑支持这一事实，但是，它被要求在无限范围内进行积分。如果将积分范围更改为函数支持的范围（即 [1,2]），则它没有问题。

Answer 2

dunif 以整数交易...

dunif(.9, 1, 2)
[1] 0

dunif(1, 1, 2)
[1] 1

dunif(1.5, 1, 2)
[1] 1

dunif(2, 1, 2)
[1] 1

dunif(2.1, 1, 2)
[1] 0

如果您只对函数使用整数然后进行插值，您会得到预期的答案。

convXY2(2)
[1] 0

convXY2(3)
[1] 1

convXY2(4)
[1] 0

z <- 3.75
remainder <- z %% 1
convXY2(floor(z))*(1-remainder) + convXY2(ceiling(z))*(remainder) 
[1] 0.25

您的函数适用于其他连续分布，因此您可能只需要在包装函数中考虑这种特定情况。或者使用 dunif 以外的其他东西，它可以根据你给出的精度有效地划分你的概率。

dunif2 <- function(x, min, max) { if (x >= min & x <= max) { return(10^-nchar(strsplit(as.character(x), "\\.")[[1]][[2]])) } else { return(0) } }

dunif2(1.45, 1, 2)
[1] 0.01

dunif2(1.4, 1, 2)
[1] 0.1