加快与天气相关的索赔的蒙特卡罗模拟速度

Question

我正在开展蒙特卡罗模拟项目，我需要帮助优化问题的某些方面。这是场景：

我们将索赔 R（与天气相关事件相关）建模为自变量，在某些假设下天气条件本身是独立的。这种情况涉及关于摩托车手的两个假设：

1. 假设 A：所有摩托车手都在同一地区行驶，受到相同的天气条件。

2. 假设B：摩托车手在不同的地区旅行，每个地区都有独立的天气条件（建模为独立的马尔可夫链）。

索赔金额（以货币单位）遵循分布 $\mathbb{P}r$，密度与 成正比： $f^*(x) = extbf{1}{\mathbb{R}} xp\left(- ta\ln( lpha + |x - x_0|) 右）$,

其中 $lpha > 0, x_0 > 0, ta > 0$ 是常数。

在 365 天的时间内，我们投保 $N$，总索赔金额 $R$ 是个人索赔的总和。

我的目标是使用蒙特卡罗方法近似条件期望 $m(s) = \mathbb{R}[R - s \mid R > s]$。

这是我编写的当前 R 代码，用于根据 $\mathbb{P}_r$ 模拟索赔金额：

rremb <- function(n, params) {
  # Extract parameters from the list
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Check parameter validity
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Invalid parameters: alpha > 0, x0 > 0, and eta > 3 are required")
  }
  
  # Normalization factor K
  K <- 2 * alpha^(1 - eta) / (eta - 1)
  
  # Inverse of F* (without normalization factor)
  F_star_inv <- function(z) {
    if (z <= 0.5) {
      return((alpha + x0 - (z * (eta - 1))^(1 / (1 - eta))))
    } else {
      return((x0 - alpha + ((1 - z) * (eta - 1))^(1 / (eta - 1))))
    }
  }
  
  # Inverse of F considering K
  F_inv <- function(y) {
    return(F_star_inv(K * y))
  }
  
  # Generate n uniform random values on (0, 1)
  u <- runif(n)
  
  # Use F_inv to generate values based on P_r
  x <- sapply(u, F_inv)
  
  return(x)
}

# Parameters
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)

# Simulate values
n <- 100000  # Adjust as needed
simulated_values <- rremb(n, params)

# Display results
print(head(simulated_values, 10))

# Optional: Plot histogram of simulated values
hist(simulated_values, breaks = 50, main = "Histogram of Simulated Values", 
     xlab = "Values", ylab = "Frequency", col = "skyblue")

虽然代码工作正常，但对于大型模拟来说，它的速度没有我想要的那么快。有办法加快速度吗？例如：

1. sapply 调用可以用更快的东西代替吗？

2. 有没有办法对F_inv函数进行向量化或优化？

谢谢您的帮助！

Answer 1

这是您的函数的两个版本，它们都快一个数量级。

• 第一个

  rremb2

  ，用

  if

  向量化

  ifelse

  指令；
• 第二个

  rremb3

  ，使用

  if

  条件上的逻辑索引进行向量化。

并且不需要两个

F_*_inv

，只需将产品

K * y

传递给

F_star_inv

即可。直接。

rremb2 <- function(n, params) {
  # Extract parameters from the list
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Check parameter validity
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Invalid parameters: alpha > 0, x0 > 0, and eta > 3 are required")
  }
  
  # Normalization factor K
  K <- 2 * alpha^(1 - eta) / (eta - 1)
  
  # Inverse of F* (without normalization factor)
  F_star_inv <- function(z) {
    ifelse(z <= 0.5,
      (alpha + x0 - (z * (eta - 1))^(1 / (1 - eta))),
      (x0 - alpha + ((1 - z) * (eta - 1))^(1 / (eta - 1)))
    )
  }
  
  # Generate n uniform random values on (0, 1)
  u <- runif(n)
  # Use F_inv to generate values based on P_r
  x <- F_star_inv(u * K)
  
  x
}
rremb3 <- function(n, params) {
  # Extract parameters from the list
  alpha <- params$alpha
  x0 <- params$x0
  eta <- params$eta
  
  # Check parameter validity
  if (alpha <= 0 | x0 <= 0 | eta <= 3) {
    stop("Invalid parameters: alpha > 0, x0 > 0, and eta > 3 are required")
  }
  
  # Normalization factor K
  K <- 2 * alpha^(1 - eta) / (eta - 1)
  
  # Inverse of F* (without normalization factor)
  F_star_inv <- function(z) {
    res <- numeric(length(z))
    i <- z <= 0.5
    res[i] <- alpha + x0 - (z[i] * (eta - 1))^(1 / (1 - eta))
    res[!i] <- x0 - alpha + ((1 - z[!i]) * (eta - 1))^(1 / (eta - 1))
    res
  }
  
  # Generate n uniform random values on (0, 1)
  u <- runif(n)
  # Use F_inv to generate values based on P_r
  F_star_inv(u * K)
}

# Parameters
params <- list(alpha = 2.0, x0 = 50.0, eta = 4.0)
# Adjust as needed
n <- 100000

# Simulate values setting the pseudo-RNG seed so that results are reproducible.
set.seed(2024)
simulated_values <- rremb(n, params)
set.seed(2024)
simulated_values2 <- rremb2(n, params)
set.seed(2024)
simulated_values3 <- rremb3(n, params)

# compare results
identical(simulated_values, simulated_values2)
#> [1] TRUE
identical(simulated_values, simulated_values3)
#> [1] TRUE

library(microbenchmark)
microbenchmark(
  original = rremb(n, params),
  rremb2 = rremb2(n, params),
  rremb3 = rremb3(n, params)
) |> print(order = "median", unit = "relative")
#> Unit: relative
#>      expr       min        lq      mean    median        uq       max neval cld
#>    rremb3  1.000000  1.000000  1.000000  1.000000  1.000000 1.0000000   100   b
#>    rremb2  1.041313  1.020075  1.017987  1.043388  1.033261 0.5482928   100   b
#>  original 11.501961 12.281259 12.611460 12.976389 13.417170 6.5469023   100  a

^{创建于 2024 年 11 月 23 日，使用 reprex v2.1.1}