r - 如何在不每次将结果行数加倍的情况下 convolve() 两个以上的分布？

Question

我正在尝试convolve()36 个 beta 发行版。结果分布是下一次连续调用的两个输入分布之一convolve()。在 every 之后convolve()，结果有 row count = (nrow(vector1)+nrow(vector2)-1)。实际上，每次调用convolve(). 这是非常低效的——它使运行时间变得非常长并且消耗大量内存。有没有办法保持行数不变？

下面的代码示例...

# Function from https://stat.ethz.ch/pipermail/r-help/2008-July/168762.html
weighted.var <- function(x, w, na.rm = FALSE) {
  if (na.rm) {
    w <- w[i <- !is.na(x)]
    x <- x[i]
  }
  sum.w <- sum(w)
  sum.w2 <- sum(w^2)
  mean.w <- sum(x * w) / sum(w)
  (sum.w / (sum.w^2 - sum.w2)) * sum(w * (x - mean.w)^2, na.rm = na.rm);
}

# Define beta distribution shape parameters.
s1a <- 3.52; s1b <- 65.35;
s2a <- 1.684; s2b <- 189.12;
s3a <- 5.696; s3b <- 32.34;
s4a <- 1.81; s4b <- 185.5;

# Define intial set of quantiles.
mQ1 <- matrix(data=seq(0,1,1/1000),ncol=1);

for (i in 1:3){
  mPDF <- matrix(data=convolve(dbeta(mQ1,s1a,s1b),rev(dbeta(mQ1,s2a,s2b)),type="open"),ncol=1L);
  print(paste(nrow(mPDF),' rows',sep=''));

  if(i < 3){
    # Calculate the merged shape parameters directly from mPDF.
    mQ2 <- matrix(data=seq(0,1L,(1L/(nrow(mPDF)-1L))),ncol=1L);
    wtMean <- weighted.mean(mQ2,mPDF);
    wtStd <- sqrt(weighted.var(mQ2,mPDF));
    s1a <- -1L * ((wtMean*(wtStd^2 + wtMean^2 - wtMean))/wtStd^2);
    s1b <- ((wtStd^2 + wtMean^2 - wtMean)*(wtMean - 1))/wtStd^2;

    s2a <- s3a; s2b <- s3b;
    mQ1 <- mQ2;
  }

} #i

标签： rconvolution