Base R函数卷积

Question

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

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

在这个 post 之后，我定义了以下函数，它将一对 pdf 作为参数，return 是它们的卷积：

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

更一般地说，这适用于许多均匀分布对的组合，但是当我尝试对一对均匀 [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

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

Answer 2

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

Base R函数卷积

Base R Function Convolution

r

convolution

probability-density