CODEWITHSUNDEEP
rrandomvectorizationrandom-seed
Danilo Perl
Danilo Perl

Reputation: 69

Differences between vectorized and non-vectorized code in R

I am new to R and am trying to vectorize an R script. However, my attempt at a vectorized version returns values ​​of mean(SAg) and sd(SAg) (see below) that are different from those returned by the original script, even if I use the same RNG (random number generator) seed.

This is the original, non-vectorized code:

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)

And this is my attempt at a vectorized implementation:

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

Can someone please point out what I am getting wrong?

Upvotes: 0

Views: 86

Answers (1)

hi_everyone
hi_everyone

Reputation: 333

(This ought to be a comment only, but I'm lacking reputation for that) I'm not entirely sure what's happening, but my guess is that the slight difference comes from the different order in which the (pseudo)random draws are used. In the vectorised code, R makes all the runif draws first and then all the rexp draws, in the non-vectorised version, it alternates between the two. Did you try to simulate it many times to rule out that one version structurally yields smaller mean/sd than the other version?

Upvotes: 2

Related Questions