Exemplos de simulação usando o método da transformação inversa:
# Limpeza das variaveis e funcoes do ambiente
rm(list=ls())
Exemplo 1: simulação de uma variável com distribuição exponencial
cdf: F(x; lambda) = 1 - exp(-lambda*x)
inversa: Q(u; lambda) = -log(1-u)/lambda
# Funcao para calcular a inversa da fda da distribuicao exponencial
# u: real; probabilidade para a qual se deseja obter o respectivo
# quantil; 0 < u < 1
# lambda: real; parametro da distribuicao: taxa de ocorrencia do evento
# por unidade de tempo, volume, etc. lambda>0
inv_fda_exp = function(u, lambda=1.0) {
return(-log(1-u)/lambda)
}
# funcao para simular m valores da distribuicao exponencial
# Argumentos:
# m: inteiro positivo indicando a quantidade de valores; m>=1
# lambda: real; parametro da distribuicao: taxa de ocorrencia do evento
# por unidade de tempo, volume, etc. lambda>0
# Saida:
# Vetor de tamanho m contendo os valores sorteados
sim_exp = function(m=1, lambda=1.0) {
if (m<1 | lambda<=0) {
stop('Erro na funcao sim_exp: verifique argumentos.')
}
# Sorteio de m valores da uniforme (0,1)
u = runif(m)
# Inversao para obter x ~ exponencial(lambda)
x = inv_fda_exp(u, lambda)
return(x)
}
m = 10000 # numero de pontos sorteados
lambda = 1 # parametro da distribuicao lambda (>0)
x = sim_exp(m,lambda)
# histograma da variaval gerada
h = hist(x, 100)
# Apenas para comparacao: pdf dos pontos medios das barras
# do histograma
fdp_mids = dexp(h$mids, lambda) # funcao pre-definida do R
plot(h$mids, h$counts/m, type='l', lwd=2)
lines(h$mids, fdp_mids/sum(fdp_mids), col='red', lwd=2)
legend(x='topright', legend=c('Fdp simulada', 'Fdp teorica'), lwd=2, col=c('black', 'red'))
Exemplo 2: simulação de uma variável com distribuição categórica (versão multidimensional da distribuição de bernoulli)
# funcao para simular m valores da distribuicao binomial
# Argumentos:
# m: inteiro positivo indicando a quantidade de vetores a
# serem gerados; m>=1
# p: vetor real indicando as probabilidades de cada elemendo x_i;
# p_i>=0, sum(p) = 1
# Saida:
# Vetor de tamanho m contendo os rótulos sorteados
sim_categ = function(m=1, p=rep(0.5,2)) {
if (m<1 | sum(p)!=1) {
stop('Erro na funcao sim_categ: verifique argumentos.')
}
fda_categ = cumsum(p) # Calculo da acumulada F
print(fda_categ) # Apenas para ilustracao - comentar esta linha
u = runif(m) # Sorteio de m valores da uniforme (0,1)
print(u[1:5]) # Apenas para ilustracao - comentar esta linha
#
# Aplicacao da inversa generalizada de F para obter x ~ categorica(p)
# Ver algoritmo nos slides do Capitulo 03 - Variaveis Aleatorias
#
# Funcao nativa do R que localiza em qual intervalo do vetor
# fda_categ cada valor do vetor u se encontra
x = findInterval(u, fda_categ)
print(x[1:5]) # Apenas para ilustracao - comentar esta linha
return(x)
}
m = 1000 # numero de pontos sorteados
p=c(0.2,0.6,0.2) # parametro da distribuicao (vetor de probabilidades das categorias)
x = sim_categ(m,p)
## [1] 0.2 0.8 1.0
## [1] 0.89078342 0.46166544 0.42575819 0.15664659 0.05746409
## [1] 2 1 1 0 0
# histograma da variaval gerada
barplot(rbind(table(x)/m, p), col=c('grey80', 'grey40'),
legend.text=c('F.p. simulada', 'F.p. teorica'), beside=TRUE)
abline(h=0)