a) Geração direta pelo método da função inversa;

Para obter uma amostra de observações geradas por uma distribuição binomial(n,p) pelo método da função inversa, é necessário primeiramente definir “faixas de de probabilidade”, ou tresholds. Tais faixas representam a probabilidade de um evento ocorrer. Suponha \(X \sim binom(n,p)\), dado um valor inteiro arbitário \(a\), uma faixa poderia ser definida como os valores contínuos compreendidos entre \(P(X=a-1)\) e \(P(X=a)\).

Em linhas práticas, gera-se um valor \(u\) de uma distribuição unif(0,1) e como critério comparativo, avalia-se em qual faixa este valor está contido. Uma forma prática de otimizar a eficiência computacional é organizar as faixas de forma cumulativa, tem-se então, para uma Binomial com n = 10 e p = 0,3:

[0-0,02] ; [0,02-0,14] ; [0,14-0,38] ; ….. ; [0,99-1]

Por exemplo, supondo que u = 0,03, o valor gerado através do método da inversa seria 1, pois caiu dentro da faixa definida por x = 1.

b) Geração por soma de Bernoullis;

Existem oturas formas para gerar uma amostra proveniente de uma distriuição binomial(n,p) através do método da transformação inversa, uma delas é utilizar o resultado de que uma distribuição binomial pode ser decomposta em uma soma de distribuição bernoullis, como segue:

\[ X_i\sim Ber(p) \implies Y = \sum_{i=1}^{n}X_i \sim binom(n,p) \]

Através do resultado acima, um caminho simples para gerar da distribuição Binomial é realizar 10 Bernoulli trials com probabilidade de sucesso p = 0,3 e somar a quantidade de sucessos, o resultado da soma é exatamente o valor gerado supondo uma distribuição Binomia(10;0,3).

Supondo que não é conhecido como realizar simulação de uma distribuição Bernoulli(0,3), pode-se gerar um valor u vindo de uma distriuição Uniforme(0,1) e comparar diretamente com a probabilidade de sucesso. Se u for menor que p = 0,3, então, o valor simulado para uma distribuição Bernoulli é de 0, caso contrário, 1.

c) Geração direta pelo rbinom;

A forma mais simples possível de simular uma distribuição Binomial(n,p) pelo R é através da função rbinom(q,n,p) onde o parâmetro q representa o tamanho da amostra simulada, n e p são análogos a definição teórica da distribuição Binomial(n,p).

Internamente, a rotina utilizada pela função é definida através da proposta feita por Kachitvichyanukul, V. and Schmeiser, B. W. (1988) no artigo Binomial random variate generation. Communications of the ACM, 31, 216–222. O autor avalia os métodos usuais na literatura para geração de variáveis aleatórias binomiais e compara em função da eficiência computacional, em termos de tempo.

Inspeção gráfica

Verificou pela Figura 1 poucas diferenças absolutas, todavia, observa-se que o padrão da Transformação Inversa e da Geração Rbinom são semelhantes, enquanto a Soma de Bernoullis apresenta leve diferença, em especial para as categorias centrais, isto é, 3, 4 e 5. Como era esperado, os valores mais observados estão em x = 3, uma indicativa da qualidade da simulação em termos de precisão.

Verificação da simulação

Verificou pela Tabela 1 pouca diferença entre a média teórica e as médias observadas para cada método, sendo esta diferença visível apenas na terceira cada decimal. Em outras palavras, os resultados encontrados sugerem uma boa precisão quanto aos métodos utilizados. Ademais, tem-se também uma proximidade numérica entre as variâncias teórica e observada, uma indicativa que não houve problema de sobredispersão. Caso contrário, se a variância fosse maior ou menor do que a teórica, seria uma evidência de possível erro na simulação.

A Tabela 2 apresenta o tempo de execução necessária para cada simulação, supondo tamanho 10.100 (dez mil e cem) e com burnin de 100 (cem) para evitar uma variação inicial indesejada, resultando em um tamanho amostral efetivo de 10.000 (dez mil). Verificou-se que as diferenças no custo computacional são baixas, indicando uma eficiência dos três metodo. Contudo, a Geração via rbinom apresentou o melhor desempenho neste quesito, como era esperado, pois considera o método mais otimal em relação ao custo computacional.

Distribuição N(0,1) como proposta

Para simular de uma distribuição Normal truncada entre -1 e 1 através de uma Normal(0,1), um caminho é implementer uma função indicadora que faça com que a distribuição normal proposta gere valores apenas entre -1 e 1, caso contrário o valor gerado recebe 0.

Para o cálculo da Taxa de Aceitação teórica, é necessário avaliar c tal que:

\[ c = \cfrac{f(x)}{g(x)} \] Onde \(f(x)\) é a função objetivo e a \(g(x)\) é a função proposta.

Pode-se então reescrever como

\[ c = \cfrac{k/\sqrt{2\pi}*exp\{-x^2/2\}*I(x)_{(-1,1)}}{1/\sqrt{2\pi}*exp\{-x^2/2\}} \\ =k*I(x)_{(-1,1)} \]

Onde k é a constante normalizadora. Sua obtenção pode ser feita através da seguinte integral:

\[ k\int_\Omega f(x) dx = 1 \implies k = (\int_{-1}^{1}f(x)dx)^{-1} \]

Sendo \(f(x)\) a densidade de uma distribuição Normal(0,1). Tem-se que

\[ k = c = 1,46 \implies T = \cfrac{1}{c} = 68\% \]

Distribuição U(-1,1) como proposta

Seguindo a mesma lógica do exercício anterior, pode-se definir: \[ c = \cfrac{k/\sqrt{2\pi}*exp\{-x^2/2\}*I(x)_{(-1,1)}}{2*I(x)_{(-1,1)}} \\ =\cfrac{2k}{\sqrt{2\pi}}*exp\{-x^2/2\}*I(x)_{(-1,1)} \\ \implies c = max\{\cfrac{2k}{\sqrt{2\pi}}*exp\{-x^2/2\}*I(x)_{(-1,1)}\} \]

Aqui neste contexto, max{argumento} denota o máximo da função argumento, dependente de \(x\). Uma maneira simples de encontrar o máximo é visualizar a função argumento, para este caso, como proporcional ao núcleo de uma normal truncada, ie:

\[ c \propto exp\{-x^2/2\}*I(x)_{(-1,1)} \]

O máximo de uma normal truncada será avaliada em sua média, ou seja, em \(x=0\), logo:

\[ c =\cfrac{2k}{\sqrt{2\pi}}*exp\{0^2/2\}*I(x)_{(-1,1)} = \cfrac{2k}{\sqrt{2\pi}} = 1,16 \\ \implies T = \cfrac{1}{c} = 85,56\% \]

Distribuição Cauchy(0,1) como proposta

Analogamente aos exercícios anteriores, pode-se escrever c tal que:

\[ c = \cfrac{k/\sqrt{2\pi}*exp\{-x^2/2\}*I(x)_{(-1,1)}}{1/\pi(1+x^2)} \\ =\cfrac{2k}{\sqrt{2\pi}}*exp\{-x^2/2\}\pi(1+x^2)*I(x)_{(-1,1)} \\ \implies c = max\{ \cfrac{2k}{\sqrt{2\pi}}*exp\{-x^2/2\}\pi(1+x^2)*I(x)_{(-1,1)} \} \]

Diferentemente das funções anteriores, o máximo da função apresentada acima não tão imediata. Posto isso, foi utilizado um recurso de otimização para encontrar o máximo global, utilizando a função optimize. COm isso, tem-se que:

\[ c = 2,22 \implies T = \cfrac{1}{c} = 44,90\% \]

Taxas de aceitação

Questão 1

set.seed(07082019)

#Definindo os parâmetros de acordo com o exercício
k = 10
p = 0.3
n = 10100

#Criando um vetor vazio para armazenar as probabilidades de sucesso
prob = vector()

#Gerando todas as possíveis probabilidades
for (i in 0:k){prob[i+1]=dbinom(i,k,p)}

#Verificando se somam 1
#sum(prob)
#De acordo

#Criando as faixas (tresholds)
faixa = cumsum(prob)

#Criando um vetor vazio para armazenar os valores gerados pela binomial
x1 = vector()

#Testando através de TRUE ou FALSE  
tic()
for (i in 1:n){
u = runif(1)
x1[i] = sum(u>faixa[1:11])
}
exectime <- toc()
tempo1 <- exectime$toc - exectime$tic

#removendo burnin
x1 = x1[101:10100]

set.seed(07082019)

# Gerar um número aleatório U.
# Determinar o menor inteiro positivo I tal que:

gera_ber = function(n,p){
u = runif(n)
x=NULL
  for (i in 1:length(u)){
    x[i] = ifelse(u[i] > p,0,1)
  }
x
}

#Definindo os parâmetros
n = 10100
p = 0.3

#Algoritmo Victor
tic()
x2 <- rep(0,n)
for (i in 1:n){
  aux <- 0
  for (j in 1:k){
    aux <- aux + gera_ber(1,p)
  }
  x2[i] <- aux
}
exectime <- toc()
tempo2 <- exectime$toc - exectime$tic

#Algoritmo Rafael
for(i in 1:n){x2[i]=sum(gera_ber(10,0.3))}

x2 = x2[101:10100]

set.seed(07082019)

#Definindo os parâmetros de acordo com o exercício
k = 10
p = 0.3
n = 10100

tic()
x3 = rbinom(n,k,p)
exectime <- toc()
tempo3 <- exectime$toc - exectime$tic

x3 = x3[101:10100]

Questão 2

set.seed(07082019)

# Obtenção do C
# Area abaixo da normal padrao [-1;1]

#Cálculo da constante normalizadora para a normal truncada
normz <- pnorm(1,0,1) - pnorm(-1,0,1)
k <- 1/normz
c <- k

#Taxa de aceitação teórica 
txteorica1 <- 1/c

#Criando a normal truncada como uma função
normtrunc <- function(x){ifelse(x<=-1|x>=1,0,dnorm(x,0,1)*k)}

#Adicionando a constante C
cdnorm <- function(x){dnorm(x)*c}

razao <- function(x){normtrunc(x)/cdnorm(x)}

#Vetor de variáveis gerada de uma normal truncada
x1 = vector()

#Valor da razão A = f(x)/c*g(x)
a1 = vector()

#valor gerado por uma uniforme para comparação
u1 = vector()


#Loop para aceitação/rejeição
i <- 1
v1 <- 1
n <- 100000
while (i<n){
v1 <- v1 + 1
u1[i] <- runif(1)
xprop <- rnorm(1)
a1[i] = normtrunc(xprop)/(dnorm(xprop)*c)
if (u1[i] <= a1[i]){x1[i]<-xprop;i<-i+1}   
}

txemp1 <- n/v1

set.seed(07082019)

# Obtenção do C
# Area abaixo da normal padrao [-1;1]

#Cálculo da constante normalizadora para a normal truncada
normz <- pnorm(1,0,1) - pnorm(-1,0,1)
k <- 1/normz
c <- 2*k/sqrt(2*pi)

#Taxa de aceitação teórica 
txteorica2 <- 1/c

#Criando a normal truncada como uma função
normtrunc <- function(x){ifelse(x<=-1|x>=1,0,dnorm(x,0,1)*k)}

#Adicionando a constante C
cdunif <- function(x){dunif(x,-1,1)*c}

razao <- function(x){normtrunc(x)/cdunif(x)}

#Vetor de variáveis gerada de uma normal truncada
x2 = vector()

#Valor da razão A = f(x)/c*g(x)
a2 = vector()

#valor gerado por uma uniforme para comparação
u2 = vector()

#Loop para aceitação/rejeição
i <- 1
v2<- 1
n <- 100000
while (i<n){
v2 <- v2 + 1
u2[i] <- runif(1)
xprop <- runif(1,-1,1)
a2[i] = normtrunc(xprop)/(cdunif(xprop))
if (u2[i] <= a2[i]){x2[i]<-xprop;i<-i+1}   
}

txemp2 <- n/v2

f = function(x){dt(x,1)}

C = function(x){return(normtrunc(x)/dt(x, 1))}


Cx = function(x){-C(x)}
cc = -optimize(Cx, c(-1, 1))$objective

f = function(x){cc*dt(x,1)}

razao = function(x){normtrunc(x)/f(x)}
      
m = 1
x3 = rep(0,n)
v3 = 0
while (m <= n) {
  x3[m] = rt(1,1)
  if(runif(1) < C(x3[m])/cc){m = m + 1}
  v3 = v3 + 1
}

txteorica3 = 1/cc
txemp3 = n/v3

Questão 3

set.seed(07082019)
c = 4/pi
tx = 1/c

i <- 1
v <- 1
n <- 100
x <- vector()
x1 <- vector()
y <- vector()
y1 <- vector()
while (i<=n){
v <- v + 1
a = runif(1,-1,1) 
b = runif(1,-1,1)
if (a^2 + b^2 <= 1){x[i]=a;y[i]=b;i<-i+1}
else {x1[i]=a;y1[i]=b}
}

txemp = n/v


plot(x1,y1, asp = 1, xlim = c(-1, 1),col="blue",ylab="y",xlab="x",main="Fig 11 - Geração de valores de um círculo centrado em zero")
points(x,y,col="red")
draw.circle(0, 0, 1, nv = 1000, col = NA, lty = 1, lwd = 1)
legend(1.2, 0.65, legend=c("Rejeitado", "Aceito"),col=c("blue", "red"), lty=1, cex=0.8,lwd=2)