R 中矢量化代码和非矢量化代码的区别

Question

我是 R 的新手，正在尝试矢量化 R 脚本。然而，我尝试的矢量化版本returns值mean(SAg)和sd(SAg)（见下文）与原始脚本返回的值不同，即使我使用相同的RNG （随机数生成器）seed.

这是原始的非矢量化代码：

set.seed(42)
#n: tamanho da carteira (quantidade de segurados
tamanho = 10000
# prob_sin é a probabilidade esperada de ocorrência de sinistros com um segurado (1_i)
prob_sin = 0.02
# "cenarios" são as quantidades de simulações de atribuição de valor para o SAg
cenarios = 100
# sev_media é o valor esperado da indenização em caso de materialização do sinistro com a i-ésima apólice
sev_media = 10000

replicacoes_ind = matrix(NA, tamanho, cenarios) # 0 e 1
replicacoes_sev = matrix(NA, tamanho, cenarios) # se a matriz anterior tiver 1, então terá entrada diferente de zero

SAg = array(NA, cenarios)

# quant_sin é um vetor que armazena a quantidade de sinistros ocorridos em cada simulação/cenário
quant_sin = array(NA, cenarios)

for (j in 1:cenarios){ ### cria todos os cenários possíveis de ocorrência
  for (i in 1:tamanho){ ### percorre toda a carteira
    u <- runif(1,0,1) # gera um número aleatório entre 0 e 1
    ifelse(u <= prob_sin, replicacoes_ind[i,j] <- 1, replicacoes_ind[i,j] <- 0)
    ifelse(u <= prob_sin, replicacoes_sev[i,j] <- rexp(1,rate=1/sev_media), replicacoes_sev[i,j] <- 0)
  }
  quant_sin[j] <- sum(replicacoes_ind[,j])
  SAg[j] <- sum(replicacoes_sev[,j])
}

mean(SAg); sd(SAg)

这是我对矢量化实现的尝试：

simulate_uniform = function(tamanho, cenarios){
  return(matrix(data=runif(tamanho*cenarios, 0, 1), nrow=tamanho, ncol=cenarios))
}

fill_sev = function(p, u, dist, media, tamanho, cenarios){
  # Matriz que armazena as indenizacoes ($) associadas a cada sinistro
  matriz = matrix(NA, tamanho, cenarios)
  
  # Indices aos quais houve sinistros
  ind <- which(u <= p, arr.ind = TRUE)
  
  # Usamos os indices acima para preencher as posicoes com a distribuicao escolhida
  if(dist == "Exponencial"){
    matriz[ind] = rexp(nrow(ind), rate = 1/media)
  } else if(dist == "Uniforme"){
    matriz[ind] = runif(nrow(ind), min = 0, max = 20000)
  } else {
    return(NULL)
  }
  
  # Para os indices restantes, substituimos por zero
  matriz[is.na(matriz)] = 0
  
  return(matriz)
}

set.seed(42)

cenarios = 100 # Numero de simulacoes
tamanho = 10000 # Carteira
prob_sin = 0.02 # Probabilidade de ocorrência de sinistros
sev_med = 10000 # Severidade média individual
distribuicao = "Exponencial" # Entre a distribuicao exponencial ou uniforme

SAg = array(NA, cenarios)

# Armazena a qtd de sinistros ocorridos em cada cenario
quant_sin = array(NA, cenarios)

u = simulate_uniform(tamanho, cenarios)

replicacoes_sev = fill_sev(prob_sin, u, distribuicao, sev_med, tamanho, cenarios)

if(is.null(replicacoes_sev)){
  stop("Somente sao permitidas as ditribuicoes exponencial e uniforme.")
} else{
  # Quantidade de sinistros em cada simulacao
  quant_sin = colSums(replicacoes_sev != 0)
  
  # Sinistro agregado em cada simulacao
  SAg = colSums(replicacoes_sev)
  
  # Media de SAg
  mean(SAg)
  
  # Desvio padrao de SAg
  sd(SAg)
}

有人可以指出我哪里错了吗？

Answer 1

（这应该只是评论，但我对此缺乏声誉）我不完全确定发生了什么，但我的猜测是细微差别来自于使用（伪）随机抽奖的不同顺序。在矢量化代码中，R 首先绘制所有 runif，然后绘制所有 rexp，在 non-vectorised 版本中，它在两者之间交替。您是否尝试多次模拟它以排除一个版本在结构上比另一个版本产生的 mean/sd 更小？

R 中矢量化代码和非矢量化代码的区别

Differences between vectorized and non-vectorized code in R

r

random

vectorization

random-seed