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

I try to compute the variance of ratio between 2 sum of Gamma distribution with shape and scale different.

Numerator :

Denominator :

To compute this variance, I am using coga library of R. Here's below the code snippet :

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,]))
}

1. The first issue is that print('var(y4) = ', var(y4[i,])) doesn't display anything and I don't understand why : below the result in R shell :

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

This first issue has been solved by just using paste0 function (see first comment below).

2. The second main issue is that I get this same following error during execution :

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

However, I can plot a density function which looks like :

What is wrong with this systematical error about multiple of replacement length ? Sorry, I begin with R language.

A first part without the grid and pdf/cdf part

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,])))
}

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

and five plots similar to this one