R: Replicating samples for a function then cumulatively adding the result

I have this code where I am using the concept of uncertainty propagation in Monte Carlo simulation where instead of supplying fixed values, I am using probability distributions as inputs.

library(triangle)

# Define the model function
model_function <- function(X, Y, Z) {
  return (X * Y * Z)
}

# Number of simulations
n <- 10000

# Generate random samples for X, Y, and Z
X_samples <- rtriangle(n, a = 7563.474, b = 9244.246, c = 8403.86)
Y_samples <- rnorm(n, mean = 35, sd = 5)
Z_samples <- rnorm (n, mean = 0.6730, sd = 0.1)

# Apply the model function to each set of inputs
direct.loss <- model_function(X_samples, Y_samples, Z_samples)

# Summary statistics
summary(direct.loss)

# Plot the distribution of outputs
hist(direct.loss, breaks = 50, col = "skyblue", main = "Histogram of direct losses", xlab = "Direct losses")

The output of this would be a probability distribution. Now, what I want to do is to have multiple runs or replications (728 times) of generating the samples for X, Y, and Z, then cumulatively adding the result of the 728 runs.

To demonstrate this, here is a rather inefficient (and possibly unrealistic setup). In this code, I have two runs. The second run being the samples generated for X.1, Y.1, and Z.1.

library(triangle)

# Define the model function
model_function <- function(X, Y, Z, X.1, Y.1, Z.1) {
  return ((X * Y * Z) + (X.1 * Y.1 * Z.1))  # The second run being added to the first run
}

# Number of simulations
n <- 10000

# Generate random samples for (X, Y, and Z) and (X.1, Y.1, and Z.1)
X_samples <- rtriangle(n, a = 7563.474, b = 9244.246, c = 8403.86)
Y_samples <- rnorm(n, mean = 35, sd = 5)
Z_samples <- rnorm (n, mean = 0.6730, sd = 0.1) 

X.1_samples <- rtriangle(n, a = 7563.474, b = 9244.246, c = 8403.86) # Note that the inputs are the same as above
Y.1_samples <- rnorm(n, mean = 35, sd = 5)
Z.1_samples <- rnorm (n, mean = 0.6730, sd = 0.1)

# Apply the model function to each set of inputs
direct.loss <- model_function(X_samples, Y_samples, Z_samples, X.1_samples, Y.1_samples, Z.1_samples)

# Summary statistics
summary(direct.loss)

# Plot the distribution of outputs
hist(direct.loss, breaks = 50, col = "skyblue", main = "Histogram of direct losses", xlab = "Direct losses")

Is there a way to do the process for 728 runs with a more efficient and simplified code?

You can use replicate() to repeat your function n times and then just add the rows.

model_function <- function() {
  n <- 10000
  X <- rtriangle(n, a = 7563.474, b = 9244.246, c = 8403.86)
  Y <- rnorm(n, mean = 35, sd = 5)
  Z <- rnorm (n, mean = 0.6730, sd = 0.1) 
  X*Y*Z
}

direct.loss <- rowSums(replicate(n=728, model_function()))

hist(direct.loss, breaks = 50, col = "skyblue", main = "Histogram of direct losses", xlab = "Direct losses")