r - r 指数分布 rv 卷积的密度函数有时会产生错误值

Question

我试图在 John Wakely 的“Coalescent Theory: An Introduction”中复制图 3.4 中的情节。它由“T total”的密度函数组成，我不会进一步定义。

我的问题是我似乎无法正确地为分布函数编写函数。

分布函数应该是 (n-1) 个随机指数分布变量的标准卷积，其中第 i 个 rv 的比率为 (i/2)。基本分布可以写成：

我想我已经正确编码了：

DistExpConv <- function(lambdas, t) {
  product = function(vector, entrance) {
    sapply(entrance, function(x){prod(vector[-x] / (vector[-x] - vector[x]))})
  }

  sapply(t, function(y){
    sum(sapply(lambdas, function(x){
      x * exp(-x * y) * product(lambdas, which(lambdas == x))
    }))
  })
}

我认为这应该给出正确的响应，至少我看不到它应该在哪里出错。然后我实现了 T total 的分布：

T_totaldist <- function(n, t) {
  lambdas = sapply(2:n, function(x){ (x-1)/2 })
  DistExpConv(lambdas , t)
}

正如所见，这个函数只是形成 n-1 个 lambda，然后将它们发送到DistExpConv.

到目前为止一切顺利。T_totaldist尝试使用曲线函数绘制 n = 100时出现我的问题：

n = c(2, 5, 10, 20, 50, 100)
col = rainbow(length(n))
for(i in 1:length(n)){
  curve(T_totaldist(n = n[i], x), 0, 14, col = col[i], add = i!=1)
}
legend(8,0.5,paste("n = ",n, sep = ""), col=col, lty=1)

由 n = 100 产生的（在本例中为粉红色）曲线跳进跳出极端负值和正值。因此，我必须得出结论，我以前的功能做错了什么。奇怪的是，另一个 n 的曲线看起来很完美，并且从 t >= 2 开始，粉红色曲线也是正确的。所以我认为我可能会为小 t 值产生一些错误？我不知道如何继续使情节变得漂亮。

我可以在曲线的 x 和 y 值中看到，例如T_totaldist(100, 0.14)输出 4.252961e+13 而T_totaldist(100, 0.28)输出 -4.982278e+10。这显然不是我想要的，因为在这些小的 t 值下密度应该非常接近于零，并且密度永远不应该是负数。

标签： r