R：复制函数的样本，然后累加结果

Question

我有这段代码，我在蒙特卡洛模拟中使用不确定性传播的概念，我使用概率分布作为输入，而不是提供固定值。

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

其输出将是概率分布。现在，我想做的是多次运行或重复（728 次）生成 X、Y 和 Z 的样本，然后累积添加 728 次运行的结果。

为了证明这一点，这是一个相当低效（并且可能不切实际的设置）。在此代码中，我运行了两次。第二次运行是为 X.1、Y.1 和 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")

有没有办法用更高效和简化的代码来完成 728 次运行的过程？

Answer 1

您可以使用

replicate()

重复您的函数

n

次，然后只需添加行。

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

---