R - Error with assignments and mapply function - same error : "number of items to replace is not a multiple of replacement length"

Question

我尝试计算形状和尺度不同的 Gamma 分布的 2 个总和之间的比率方差。

分子：

分母：

为了计算这个方差，我使用了 R 的 coga 库。下面是代码片段：

library(coga)

options(max.print=100000000)
my_data <- read.delim("Array_total.txt", header = FALSE, sep = " ")
array_2D <- array(my_data)

z_ph <- c(0.9595, 1.087, 1.2395, 1.45, 1.688)
b_sp <- c(1.42904922, 1.52601862, 1.63866958, 1.78259615, 1.91956918)
b_ph <- c(sqrt(1+z_ph))

y0 <- array(1000, dim=c(5,36))
y1 <- array(0, dim=c(5,36))
y2 <- array(0, dim=c(5,36))

y2_sp <- array(0, dim=c(5,36))
y2_ph <- array(0, dim=c(5,36))

y3 <- array(0, dim=c(5,36))
y4 <- array(0, dim=c(5,36))

y1 <- ((array_2D[,1])+1)/2
y2 <- 2*array_2D[,2]

for (i in 1:5) {
y2_sp[i] <- 2*array_2D[,2]*(b_sp[i]/b_ph[i])
y2_ph[i] <- 2*array_2D[,2]
}

for (i in 1:5) {
  y3[i] <- mapply(rcoga, y0[i], y1[i], y2_sp[i]) / mapply(rcoga, y0[i], y1[i], y2_ph[i])
  y4[i] <- as.double(unlist(y3[i]))
}

# do grid
grid <- seq(0, 15, length.out=36000)
## calculate pdf and cdf
pdf <- array(0, dim=c(5))
for (i in 1:5) {
pdf[i] <- mapply(dcoga, grid, shape=y1[i], rate= y2_sp[i]) / mapply(dcoga, grid, shape=y1[i], rate= y2_ph[i])
}
## plot pdf
plot(density(y4[1,]), col="blue")
for (i in 1:5) {
  print('var(y4) = ', var(y4[i,]))
}

第一个问题是 print('var(y4) = ', var(y4[i,])) 没有显示任何内容，我不明白为什么：在 R shell 中的结果下方：

> source("exact_example.R")
[1] "var(y4) = "
[1] "var(y4) = "
[1] "var(y4) = "
[1] "var(y4) = "
[1] "var(y4) = "

第一个问题已通过使用 paste0 函数解决（见下面的第一条评论）。

第二个主要问题是我在执行过程中遇到了同样的以下错误：

Warning messages:
1: In y2_sp[i] <- 2 * array_2D[, 2] * (b_sp[i]/b_ph[i]) :
   number of items to replace is not a multiple of replacement length
2: In y2_ph[i] <- 2 * array_2D[, 2] :
   number of items to replace is not a multiple of replacement length

Warning messages:
1: In y3[i] <- mapply(rcoga, y0[i], y1[i], y2_sp[i])/mapply(rcoga,  :
      number of items to replace is not a multiple of replacement length
2: In y3[i] <- mapply(rcoga, y0[i], y1[i], y2_sp[i])/mapply(rcoga,  :
   number of items to replace is not a multiple of replacement length

Warning messages:
1: In pdf[i] <- mapply(dcoga, grid, shape = y1[i], rate = y2_sp[i])/mapply(dcoga,  :
   number of items to replace is not a multiple of replacement length
2: In cdf[i] <- mapply(pcoga, grid, shape = y1[i], rate = y2_sp[i])/mapply(pcoga,  :
   number of items to replace is not a multiple of replacement length

但是，我可以绘制如下所示的密度函数：

这个关于 multiple of replacement length 的系统错误有什么问题？抱歉，我从 R 语言开始。

Answer 1

没有 grid 和 pdf/cdf 部分的第一部分

library(coga)
library(dplyr)

# options(max.print = 100000000)

# read data ----
my_data <- read.delim("Array_total.txt", header = FALSE, sep = " ")

data <- my_data %>% 
  transmute(shape = V1,
            across(2:6, ~ .x, .names = "rate_{.col}"))


# set parameters ----
z_ph <- c(0.9595, 1.087, 1.2395, 1.45, 1.688)
b_sp <- c(1.42904922, 1.52601862, 1.63866958, 1.78259615, 1.91956918)
b_ph <- sqrt(1+z_ph)

y2_sp <- array(0, dim=c(5,36))
y2_ph <- array(0, dim=c(5,36))

# calculate l ----
y1 <- (data$shape + 1) / 2

# prepare output array ----
y4 <- matrix(0, nrow = 5, ncol = 1000)


# loop over different rate buckets ----
for (i in 1:5) {
  y2_sp[i,] <- 2 * data[, i + 1] * (b_sp[i] / b_ph[i])^2
  
  y2_ph[i,] <- 2 * data[, i + 1]
  
  # y3 <- rcoga(1000, y1, y2_sp[i,]) / rcoga(1000, y1, y2_ph[i,])
  
  y4[i,] <- rcoga(1000, y1, y2_sp[i,]) / rcoga(1000, y1, y2_ph[i,])
  
  # plot graph ----
  plot(density(y4[i,]), col="blue")
  
  # print variance ----
  print(paste0('var(y4) = ', var(y4[i,])))
}

这个returns

[1] "var(y4) = 6.71357918252685e-05"
[1] "var(y4) = 6.64394691819802e-05"
[1] "var(y4) = 5.82709872201527e-05"
[1] "var(y4) = 5.13592254057074e-05"
[1] "var(y4) = 4.43354125364343e-05"

还有五个与此相似的地块

R - Error with assignments and mapply function - same error : "number of items to replace is not a multiple of replacement length"

R - Error with assignments and mapply function - same error : "number of items to replace is not a multiple of replacement length"

r

distribution

gamma-distribution

variance