Lista do livro simulation do Sheldon Ross

1 Capítulo 3

O Capítulo 3 do livro simulation do Sheldon Ross fala sobre a geração de números pseudo-aleatórios e sobre o Método de Monte Carlo para aproximar a resolução de integrais.

1.1 Questões

1. Se \(x_0 = 5\) e \[x_n = 3x_{n-1} \ mod \ 150\] encontre \(x_1, ..., x_{10}.\)

2. Se \(x_0 = 3\) e \[x_n = (5x_{n-1} + 7) \ mod \ 200\] encontre \(x_1, ..., x_{10}.\)

Nos exercícios 3-9 use simulação para aproximar as integrais. Compare suas estimativas com o valor exato, se for conhecido.

3. \(\int_0^1 exp(e^x) dx\)

4. \(\int_0^1 (1-x^2)^{3/2} dx\)

5. \(\int_{-2}^{2} e^{x+x^2} dx\)

6. \(\int_{0}^{\infty} x(1 + x^2)^{-2} dx\)

7. \(\int_{-\infty}^{\infty} e^{-x^2} dx\)

8. \(\int_{0}^{1}\int_{0}^{1} e^{(x+y)^2} dy dx\)

1.2 Respostas

1

f1 <- function(n){
  x <- numeric(n+1) # Vetor para armazenar os valores
  x[1] <- 5 # semente
  
  # loop para gerar os números pseudo-aleatórios
  for (i in 2:(n+1)){
    x[i] <- (3*x[i-1]) %% 150
  }
  
  # retornará os números gerados sem a semente
  return(x[2:(n+1)])
}

# 10 primeiros números pseudo-aleatórios solicitados
f1(10)

##  [1]  15  45 135 105  15  45 135 105  15  45

Observe que o período desse gerador é muito curto, havendo repetição de valores rapidamente. Para que os números gerados tenham boas propriedades, é necessário escolher cuidadosamente os valores para os parâmetros.

2

f2 <- function(n){
  x <- numeric(n+1) # Vetor para armazenar os valores
  x[1] <- 3 # semente
  
  # loop para gerar os números pseudo-aleatórios
  for (i in 2:(n+1)){
    x[i] <- ((5*x[i-1])+7) %% 200
  }
  
  # retornará os números gerados sem a semente
  return(x[2:(n+1)])
}

# 10 primeiros números pseudo-aleatórios solicitados
f2(10)

##  [1]  22 117 192 167  42  17  92  67 142 117

Observe que o período desse gerador é muito curto, havendo repetição de valores rapidamente. Para que os números gerados tenham boas propriedades, é necessário escolher cuidadosamente os valores para os parâmetros.

3

# Definindo a função que será integrada
h3 <- function(x){
  exp(exp(x))
}

# Método de Monte Carlo
f3 <- function(funcao=h, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f3(h3); I

## [1] 6.326876

# Valor exato da integral solicitada
Iexato <- integrate(h3, 0, 1)$value; Iexato

## [1] 6.316564

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.01031172

4

# Definindo a função que será integrada
h4 <- function(x){
  (1 - x^2)^(3/2)
}

# Método de Monte Carlo
f4 <- function(funcao=h, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f4(h4); I

## [1] 0.5882481

# Valor exato da integral solicitada
Iexato <- integrate(h4, 0, 1)$value; Iexato

## [1] 0.5890486

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.0008005536

5

# Definindo a função que será integrada
h5 <- function(x){
  exp(x+x^2)
}

# Fazendo y = (x-a)/(b-a), dy = dx/(b-a) para colocar no intervalo (0,1)
g5 <- function(y, a, b){
  h5(a+(b-a)*y)*(b-a)
}

# Método de Monte Carlo
f5 <- function(funcao, a, b, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U, a, b) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f5(g5, -2, 2); I

## [1] 92.31046

# Valor exato da integral solicitada
Iexato <- integrate(h5, -2, 2)$value; Iexato

## [1] 93.16275

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.852291

6

# Definindo a função que será integrada
h6 <- function(x){
  x*(1+x^2)^(-2)
}

# Fazendo y = 1/(x+1), dy = -dx/(x+1)^2 para colocar no intervalo (0,1)
g6 <- function(y){
  h6((1/y)-1)/(y^2)
}

# Método de Monte Carlo
f6 <- function(funcao, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f6(g6); I

## [1] 0.4976929

# Valor exato da integral solicitada
Iexato <- integrate(h6, 0, Inf)$value; Iexato

## [1] 0.5

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.002307126

7

# Definindo a função que será integrada
h7 <- function(x){
  exp(-(x)^2)
}

# Fazendo y = 1/(1+exp(-x)), dy = exp(-x)/(1+exp(-x))^2 para colocar no intervalo (0,1)
g7 <- function(y){
  h7(log(y/(1-y)))/(y-y^2)
}

# Método de Monte Carlo
f7 <- function(funcao, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f7(g7); I

## [1] 1.767267

# Valor exato da integral solicitada
Iexato <- integrate(h7, -Inf, Inf)$value; Iexato

## [1] 1.772454

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.005186576

8

# Definindo a função que será integrada
h8 <- function(x, y){
  exp((x+y)^2)
}

# Método de Monte Carlo
f8 <- function(funcao, n=100000){
  U1 <- runif(n) # Gerando números da distribuição U(0,1)
  U2 <- runif(n)
  H <- funcao(U1, U2) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f8(h8); I

## [1] 4.902097

# Valor exato da integral solicitada
library(cubature)

h8_1 <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  
  exp((x1+x2)^2)
}

Iexato <- adaptIntegrate(h8_1, c(0, 0), c(1, 1))$integral; Iexato

## [1] 4.899159

# Erro da aproximação feita
abs(I - Iexato)

## [1] 0.002938178

2 Capítulo 4

O capítulo 4 do livro simulation do Sheldon Ross fala sobre a geração de variáveis aleatória discretas por meio de métodos como a tranformada inversa, aceitação e rejeição, entre outros.

2.1 Questões

1. Escreva um programa que gere n valores da distribuição de massa de probabilidade p1 = \(\frac{1}{3}\) e p2 = \(\frac{2}{3}\).

   1. Considere \(n = 100\) e execute o programa e detemine a proporção de valores iguais a \(1\);
   2. Repita (a) para \(n = 1000\).
   3. Repita (a) para \(n = 10000\).

2. Escreva um programa que, dada uma função de massa de probabilidade \(\{pj , j = 1, . . . , n\}\) como entrada, retorne como uma saída os valores da variável aleatória com a função de massa de probabilidade informada na entrada.

3. Apresente um algoritmo eficiente para simular valores de uma variável aleatória X tal que \[P\{X = 1\} = 0, 3,\, P\{X = 2\} = 0, 2, \, P\{X = 3\} = 0, 35,\, P\{X = 4\} = 0, 15.\]

4. Use um procedimento eficiente, juntamente com o the text’s random number sequence (cap 3, exercício 14), para gerar uma sequencia de 25 variáveis aleatórias Bernoulli independentes, cada uma com parâmetro \(p = 0, 8\). Quantos números aleatórios foram necessários.

5. A função de massa de probabilidade da distribuição binomial negativa com parâmetros \((r, p)\), em quer é um inteiro positivo e \(0 < p < 1\) é dada por \[\frac{(j-1)!}{(j-r)!(r-1)!}p^r(1-p)^{j-r}, \, j=r,r+1, ...\]

   1. Use a relação entre as variáveis aleatórias binomial negativa e geométrica e o resultado do Exemplo 4d do livro texto para obter um algoritmo para simular desta distribuição.
   2. Verifique a relação \[p_{j+1}=\frac{j(1-p)}{j+1-r}p_j\]
   3. Use a relação apresetada em (b) para obter um segundo algoritmo para gerar variáveis aleatórias binomial negativa.
   4. Use a interpretação da distribuição binomial negativa como o número de tentativas até a obtenção deum total de r sucessos quando os resultados de cada tentativa são indenpendentes com probabilidade de sucesso p, para obter ainda uma outra abordagem para gerar valores desta variável aleatória.

6. Apresente dois métodos para gerar uma variável aleatória X tal que \[P\{X = i\} = \frac{e^{-\lambda}\lambda^i/i!}{\sum_{j = 0}^ke^{-\lambda}\lambda^j/j!}, i = 0, ..., k\]

7. Se \(Z\) é uma variável aleatória normal padrão, mostre que \[E[|Z|]=\bigg(\frac{2}{\pi}\bigg)^\frac{1}{2}=0,798\]

8. Um par de dados são continuamente lançados até que todos os possíveis resultados; \(2, 3, 4,...,12\) tenham ocorrido pelo menos uma vez. Desenvolva um estudo de simulação para estimar o número esperado de números de lançamentos de dados necessários.

9. Suponha que a variável aleatória X pode tomar qualquer um dos valores \(1, . . . , 10\) com respectivas probabilidades \(0,06, 0,06, 0,06, 0,06, 0,06, 0,15, 0,13, 0,14, 0,15, 0,13\). Use a abordagem da composição para apresentar um algoritmo que gere ros valores de X. Use o text’s random number sequence para gerar X.

10. Suponha que \(0 ≤ λn ≤ λ\) para todo \(n ≥ 1\). Considere o seguinte algoritmo para gerar uma variável aleatória tendo taxa de hazard discreta \(\{λ\}\).
    PASSO 1: S = 0.
    PASSO 2: Gere U e faça \(Y = Int\bigg(\frac{log(U)}{log(1-\lambda)}\bigg)+1\)
    PASSO 3: \(S = S + Y\).
    PASSO 4: Gere U.
    PASSO 5: Se \(U\leq \lambda_s/ \lambda\), faça \(X = Y\) e pare. Caso contrário, vá para 2.

    1. Qual a distribuição de Y no Passo 2?
    2. Explique o que o algortimo está fazendo.
    3. Argumente sobre X ser uma variável aleatória com taxa de hazard discreta \(\{λ\}\).

2.2 Respostas

1

Letra (a)

f1 <- function(n, valor = c(1, 0)){
  x <- numeric(n) # vetor que vai armazenar o números
  
  # loop que vai gerar os números
  for (i in (1:n)){
    u <- runif(1)
    if (u < 1/3) x[i] <- valor[1]
    else x[i] <- valor[2]
  }
  
  return(x) # retorna o vetor x
}


# Proporção dos números gerados
n <- 100
prop <- table(f1(n))/n; prop

## 
##    0    1 
## 0.73 0.27

Letra (b)

# Proporção dos números gerados
n <- 1000
prop <- table(f1(n))/n; prop

## 
##     0     1 
## 0.682 0.318

Letra (c)

# Proporção dos números gerados
n <- 10000
prop <- table(f1(n))/n; prop

## 
##      0      1 
## 0.6676 0.3324

2

A função que criei sempre gera números de 1 até o tamanho do vetor de probabilidades. Logo, \[P(X = 1) = p_1, P(X = 2) = p_2, ..., P(X = n) = p_n \]

f2 <- function(n, probabilidades){
  v <- 1:length(probabilidades) # Valores que a VA X toma
  u <- runif(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  
  for (j in 1:n){
    i <- 1
    prob <- probabilidades[i] # probabilidade de posição i
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[j] < F){
        X[j] <- v[i]
        pare = "OK"
      }
      
      prob <- probabilidades[i+1]
      F <- F + prob
      i <- i + 1
    }
  }
  return(X)
}

n <- 1000
p <- c(0.3, 0.1, 0.2, 0.4) # A soma deve dar 1
x <- f2(n, p); x

##    [1] 3 2 3 1 4 2 1 4 1 2 4 3 4 4 4 1 1 2 3 1 4 1 4 1 4 1 4 4 4 4 1 4 1 2 1 4 1
##   [38] 1 4 2 1 3 3 1 4 3 1 1 3 1 1 1 3 4 3 1 1 1 4 4 1 2 3 4 4 4 4 1 2 3 2 4 2 1
##   [75] 3 3 3 3 1 4 4 4 4 2 4 4 4 4 4 1 4 4 4 3 4 4 1 3 2 4 3 4 2 1 3 3 4 4 3 4 1
##  [112] 4 2 4 2 1 2 4 1 3 3 4 4 1 1 4 3 4 3 1 1 1 3 1 4 1 1 3 2 4 2 4 1 3 2 1 4 1
##  [149] 4 1 3 1 3 3 4 1 4 4 4 4 1 3 3 1 3 1 3 4 3 4 1 4 1 3 4 1 4 4 1 4 3 1 4 1 1
##  [186] 3 1 4 1 1 3 2 3 1 2 3 1 4 4 4 1 4 4 3 4 3 2 1 3 3 4 1 1 4 3 4 2 1 1 2 2 3
##  [223] 1 3 2 4 2 1 2 1 3 3 1 4 4 1 1 1 4 4 1 2 1 4 1 1 4 4 4 4 4 4 4 1 3 4 4 2 4
##  [260] 4 4 4 4 1 4 1 4 4 4 2 2 3 4 4 4 2 4 3 3 2 4 3 1 3 1 1 2 4 1 4 1 3 2 3 3 3
##  [297] 4 3 4 4 4 4 4 4 3 3 3 4 4 2 1 3 1 3 4 1 4 3 4 1 3 3 4 1 3 1 1 4 1 4 3 4 1
##  [334] 4 2 3 4 1 1 4 2 2 4 4 1 3 1 1 2 4 3 4 3 4 4 1 3 1 3 1 3 1 1 4 4 1 4 1 4 1
##  [371] 1 3 4 1 4 4 4 3 4 4 4 1 2 3 4 1 4 4 4 4 3 1 4 3 1 3 2 1 4 3 1 1 2 1 3 1 1
##  [408] 1 2 3 1 1 1 3 4 4 4 4 3 1 1 4 4 2 3 1 2 1 3 4 3 4 4 4 4 3 4 4 1 3 4 3 1 1
##  [445] 1 1 2 1 4 2 1 2 4 4 3 4 4 3 1 3 4 3 1 2 4 4 4 1 1 3 4 1 1 1 4 2 3 4 4 4 1
##  [482] 2 4 1 4 4 2 4 1 4 4 1 2 2 2 2 3 4 4 1 4 4 3 1 1 1 1 1 4 4 4 1 1 4 3 1 4 4
##  [519] 4 1 3 3 2 1 1 1 1 1 3 4 3 4 4 3 4 1 4 3 2 1 4 1 1 1 1 1 1 1 1 4 3 1 1 1 1
##  [556] 1 4 4 1 3 1 4 2 2 4 1 1 4 4 4 4 4 1 1 2 4 4 4 4 4 4 3 1 1 2 3 4 1 2 4 4 4
##  [593] 4 3 4 1 1 1 1 4 4 4 1 4 4 1 1 1 4 3 2 2 4 4 1 3 4 1 3 4 2 1 1 4 3 1 4 3 4
##  [630] 1 4 3 1 3 4 4 4 4 4 1 1 3 3 3 1 1 4 3 1 4 1 4 2 1 4 4 3 2 3 1 1 1 1 1 4 1
##  [667] 3 4 1 3 1 3 4 4 3 4 3 2 4 2 1 4 2 1 4 4 4 3 3 2 1 2 1 2 1 4 4 1 1 1 3 1 1
##  [704] 1 4 4 4 4 3 3 3 4 1 3 4 2 3 3 4 3 4 4 1 1 3 3 4 4 4 4 1 1 4 4 2 4 3 1 4 4
##  [741] 4 3 4 3 2 1 1 4 4 1 3 4 1 4 4 4 3 4 4 4 4 2 4 2 1 1 2 1 3 1 3 4 1 1 3 1 1
##  [778] 1 1 1 3 4 4 3 1 4 2 1 4 3 3 4 1 4 4 4 4 2 1 4 3 4 1 4 1 4 4 3 1 1 4 3 1 1
##  [815] 4 3 4 3 4 4 1 3 1 1 3 1 3 1 1 4 2 1 1 1 4 1 2 3 4 3 1 4 3 4 4 4 1 4 1 1 1
##  [852] 1 4 4 3 4 4 4 4 3 1 1 2 4 4 3 1 4 1 4 3 4 4 1 3 3 1 3 3 1 2 2 3 3 3 4 3 1
##  [889] 3 4 4 4 2 3 4 1 4 4 2 3 2 4 4 3 1 4 1 1 1 1 4 1 1 1 3 3 1 1 4 4 1 1 1 4 4
##  [926] 1 2 3 1 4 4 4 1 1 1 4 4 1 4 1 4 1 2 3 3 4 4 2 2 4 1 4 4 3 4 4 3 1 4 4 4 4
##  [963] 1 2 1 1 4 3 1 4 4 1 3 3 1 4 2 4 3 2 1 4 3 3 4 1 1 3 1 4 1 4 1 4 3 1 1 3 1
## [1000] 4

table(x)/n

## x
##     1     2     3     4 
## 0.322 0.099 0.202 0.377

3

Esse gerador é mais eficiente por gerar primeiro os números com maior probabilidade.

f3 <- function(n){
  x <- numeric(n) # vetor que vai armazenar o números
  
  # loop que vai gerar os números
  for (i in (1:n)){
    u <- runif(1)
    if (u < 0.35) x[i] <- 3
    else {
      if (u < 0.35 + 0.3) x[i] <- 1
      else{
        if (u < 0.35 + 0.3 + 0.2) x[i] <- 2
        else{x[i] <- 4}
      }
    }
  }
  
  return(x) # retorna o vetor x
}

# Proporção dos valores
n <- 10000
table(f3(n))/n

## 
##      1      2      3      4 
## 0.2989 0.1971 0.3532 0.1508

4

O método the text’s random number sequence é dado pelo seguinte código.

trns <- function(n){
  x <- numeric(n+2) # Vetor para armazenar os valores
  x[1] <- 23 # semente
  x[2] <- 66
  
  # loop para gerar os números pseudo-aleatórios
  for (i in 3:(n+2)){
    x[i] <- (3*x[i-1] + 5*x[i-2]) %% 100
  }
  
  un <- x[3:(n+2)]/100
  
  # retornará os números gerados sem a semente
  return(un)
}

Vamos gerar uma binomial(m, p) pelo método da transformada inversa utilizando the text’s random number sequence para gerar as uniformes no intervalo de 0 até 1.

f4 <- function(n, m, p){
  u <- trns(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  
  for (i in 1:n){
    a <- 0 
    prob <- (1-p)^m # probabilidade quando x = 0
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[i] < F){
        X[i] <- a
        pare = "OK"
      }
      
      a <- a + 1
      prob <- ((m-a+1)/a)*((prob*p)/(1-p))
      F <- F + prob
    }
  }
  return(X)
}

n <- 100
m <- 10
p <- 0.4
npa <- f4(n, m, p)
hist(npa, main = "Binomial")

media <- mean(npa)
mu <- m*p
media; mu

## [1] 3.91

## [1] 4

variancia <- var(npa)
sigma2 <- m*p*(1-p)
variancia; sigma2

## [1] 2.062525

## [1] 2.4

5

Letra (a)

O exemplo 4d do livro fala que a variável aleatória \[X = Int\bigg(\frac{log(U)}{log(q)}\bigg) + 1\] é uma geométrica de parâmetro \(p\). Sejam \(X_1, ..., X_r\) variáveis aleatórias i.i.d. com \(X_i \thicksim G_0(p)\). Desse modo, \(\sum_{i=1}^rX_i \thicksim BN(r, p)\).

f5a <- function(n, r, p){
  X <- numeric(n) # Vetor para armazenar os valores da VA
  
  for(i in 1:n){
    u <- runif(r)
    x <- trunc(log(u)/log(1 - p)) + 1 # x armazena r valores da Geométrica(p)
    y <- sum(x) # y é um valor da BinomialNegativa(r, p)
    
    X[i] <- y
  }
  return(X)
}

n <- 100
r <- 5
p <- 0.3

npa <- f5a(n, r, p)
hist(npa, main = "Binomial Negativa")

mu <- r/p
media <- mean(npa)
mu; media

## [1] 16.66667

## [1] 15.15

sigma2 <- (r*(1-p))/(p^2)
variancia <- var(npa)
sigma2; variancia

## [1] 38.88889

## [1] 23.78535

Letra (b)

n <- 1000
r <- 5
p <- 0.3
j <- r # quando j é igual a r (número de tentativas igual o número de sucessos)

valores <- f5a(n, r, p)
tab <- table(valores)/1000 # proporção do número de tentativas
names(tab) <- NULL

pj <- tab[1] # valor do j quando o número de tentativas foi igual o número de sucessos

pjmais1 <- ((j*(1-p))/(j+1-r))*pj # probabilidade do número de tentativas ser igual a r+1 (houve um único fracasso).

pjmais1; tab[2]

## [1] 0.0105

## [1] 0.007

erro <- abs(pjmais1 - tab[2]); erro

## [1] 0.0035

Letra (c)

f5c <- function(n, r, p){
  u <- runif(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  
  for (i in 1:n){
    j <- r # primeiro caso em que o número de tentativas é igual o número de sucessos
    prob <- p^r # probabilidade quando j = r
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[i] < F){
        X[i] <- j
        pare = "OK"
      }
      
      prob <- ((j*(1-p))/(j+1-r))*prob
      F <- F + prob
      j <- j + 1
    }
  }
  return(X)
}

n <- 100
r <- 6
p <- 0.5
npa <- f5c(n, r, p)
hist(npa, main = "Binomial Negativa")

media <- mean(npa)
mu <- r/p
media; mu

## [1] 11.74

## [1] 12

variancia <- var(npa)
sigma2 <- (r*(1-p))/(p^2)
variancia; sigma2

## [1] 11.28525

## [1] 12

Letra (d)

f5d <- function(n, r, p){
  u <- runif(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  
  for (i in 1:n){
    a <- 0 
    prob <- p^r # probabilidade quando x = 0
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[i] < F){
        X[i] <- a
        pare = "OK"
      }
      
      prob <- prob*(1-p)*((r+a)/(a+1))
      F <- F + prob
      a <- a + 1
    }
  }
  return(X)
}

n <- 100
r <- 6
p <- 0.4
npa <- f5d(n, r, p)
hist(npa, main = "Binomial Negativa")

media <- mean(npa)
mu <- (r*(1-p))/p
media; mu

## [1] 8.48

## [1] 9

variancia <- var(npa)
sigma2 <- (r*(1-p))/(p^2)
variancia; sigma2

## [1] 19.88848

## [1] 22.5

6

Transformada Inversa

f6a <- function(n, lambda, k){
  u <- runif(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  
  termos <- numeric(k)
  for (j in 0:k){
    termos[j] <- exp(-lambda)*(lambda^j)/factorial(j)
  }
  soma <- sum(termos)
  
  for (i in 1:n){
    a <- 0 
    prob <- exp(-lambda)/soma # probabilidade quando x = 0
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[i] < F){
        X[i] <- a
        pare = "OK"
      }
      
      a <- a + 1
      prob <- prob*lambda/a
      F <- F + prob
      
      if (a > k) break
    }
  }
  return(X)
}

n <- 100
lambda <- 5
k <- 10
npa <- f6a(n, lambda, k)
hist(npa)

media <- mean(npa)
variancia <- var(npa)
media; variancia

## [1] 4.84

## [1] 4.479192

Aceitação e Rejeição

# Valores exatos da fp
f6b <- function(x, lambda, k){
  termos <- numeric(k+1)
  for (j in 0:k){
    termos[j+1] <- exp(-lambda)*(lambda^j)/factorial(j)
  }
  soma <- sum(termos)
  
  y <- (exp(-lambda)*(lambda^x)/factorial(x))/soma
  return(y)
}

# Função que faz o método da aceitação e rejeição
g6b <- function(n, lambda, k){
  i <- 0:k # pontos em que a fp será calculada
  pY <- f6b(i, lambda, k) # a fp calculada nos pontos
  qY <- 1/length(pY) # uniforme discreta
  c <- max(pY/qY) # valores de c que deixam qY acima de pY
  
  x <- numeric(n) # vai receber os valores que a VA toma
  i <- 1
  
  while (i <= n){ # método aceitação e rejeição
    y <- trunc(length(pY)*runif(1)) + 1
    u <- runif(1)
    
    alpha <- pY[y] / (c*qY)
    
    if (u < alpha){
      x[i] <- y
      i <- i + 1
    }
  }
  return(x)
}

n <- 100
lambda <- 5
k <- 10
npa <- g6b(n, lambda, k)
hist(npa)

media <- mean(npa)
variancia <- var(npa)
media; variancia

## [1] 5.65

## [1] 4.411616

7

# Definindo a função que será integrada
h7 <- function(x){
  abs(x)*(1/sqrt(2*pi))*exp((-(x^2))/2)
}

# Fazendo y = 1/(1+exp(-x)), dy = exp(-x)/(1+exp(-x))^2 para colocar no intervalo (0,1)
g7 <- function(y){
  h7(log(y/(1-y)))/(y-y^2)
}

# Método de Monte Carlo
f7 <- function(funcao, n=100000){
  U <- runif(n) # Gerando números da distribuição U(0,1)
  H <- funcao(U) # Cálculo da função nos números gerados
  
  return(mean(H)) # Cálculo da média, que resulta na integral solicitada
}

I <- f7(g7)
I; (2/pi)^(1/2)

## [1] 0.798064

## [1] 0.7978846

8

f8 <- function(n){
  v <- 2:12 # Valores que a VA X toma
  u <- runif(n) # números aleatórios entre 0 e 1
  X <- numeric(n) # Vai armazenar os valores da VA
  probabilidades <- (1/36)*c(1, 2, 3, 4, 5, 6, 5, 4, 3, 2, 1)
  
  j <- 1
  while (!all(v %in% X)){
    i <- 1
    prob <- probabilidades[i] # probabilidade de posição i
    F <- prob # probabilidade acumulada
    
    pare = "NO"
    while(pare == "NO"){ # loop para gerar os valores de X pelo método da transformada inversa
      if(u[j] < F){
        X[j] <- v[i]
        pare = "OK"
      }
      
      prob <- probabilidades[i+1]
      F <- F + prob
      i <- i + 1
    }
    j <- j + 1
  }
  x_limp <- X[X != 0]
  return(x_limp)
}

n <- 100
medias <- numeric(n)

for (i in 1:n){
  npa <- f8(1000) # esse 1000 é um teto para o número de elementos de cada amostra
  medias[i] <- mean(npa) # media de cada amostra
}

mean(medias) # média das médias

## [1] 7.004984

9

f <- function(n){
  x <- numeric(n)
  
  for (i in 1:n){
    u1 <- runif(1)
    u2 <- runif(1)
    
    if (u1 < 0.1423){
      x[i] <- sample(c(7, 10), 1)
      }
    else {
      if (u1 < 0.3){
        x[i] <- sample(c(7, 8, 10), 1)
      }
      else{
        if (u1 < 0.5){
          x[i] <- trunc(5*u2) + 6
        }
        else{
          x[i] <- trunc(10*u2) + 1
        }
      }
    }
  }
  return(x)
}

n <- 1000
npa <- f(n)
table(npa)/n

## npa
##     1     2     3     4     5     6     7     8     9    10 
## 0.044 0.058 0.053 0.045 0.051 0.095 0.181 0.140 0.094 0.239

npa

##    [1]  7  4  9  7  8 10  9  6  7  4  3 10  9  7  1 10  8  8  8  8  8  8  7  9
##   [25]  9  5  8  7  7  3 10 10  3  1  6  8  9  5  1  8 10  4 10  8  7 10  7  2
##   [49] 10  7  1 10  8  6  3 10 10 10  7  7 10 10  7  7  1  6  9 10  5 10  7 10
##   [73]  3  1  5  1  9 10  1  6  9  5  5  7 10  6  8  9  1  7  2  7  5  6  5  5
##   [97]  7  7  9 10  8 10  5  7  9 10 10  8  9  7 10  6  9  7 10 10  4  8  5 10
##  [121]  7  6  8  1  2  3  8 10  8  4  8 10  5  6  6  7  2 10  8 10  6  8  1  9
##  [145] 10  3  7  3  5  3  7  6 10  7  1 10  6  6  8 10  8 10  6  9  4  2  6  9
##  [169] 10  7  5  3  3 10  6 10 10 10 10  9  3  7  8  7 10  8  6  6  2 10 10  8
##  [193]  4  6  1  9  1  1  8 10  8  8  9  7  6 10  6  3  7  9 10 10  5  8  6  8
##  [217] 10  4 10  7  5 10 10  7  7  7  7  7  8  2  7 10  8  7  2  5  3  8  5  9
##  [241] 10 10  8 10 10  6  2  7  8  6  2  1 10  7  8  1 10  9  7  9  7  8  5  9
##  [265] 10  9  7  9 10  8  1  9  3  7  7  4  7 10 10  7 10  7  3  7  7  8 10  7
##  [289]  8  8  8 10  8  6  1  7 10  7 10  6  2  4  5  7  5 10  7  7 10 10  2  1
##  [313]  7  9  4  4  6  2 10  4  9  3  3 10  3  6  8  8  7  7  7  2 10  6  6  7
##  [337]  2 10  8 10  5  6  5  6  8  9  3 10  6 10 10 10  6  3  4 10 10  7 10 10
##  [361] 10  2  1  2  3  6 10  7  6  2  9  6  7  7  9 10  3  9  9  7 10  3  5 10
##  [385]  5  9  4  8  2  7 10  8  4  7 10 10  7  4  8 10 10  6  7  8  7 10  9  9
##  [409]  8  6  9  7  3 10  7  3  9 10  5  6  3  8  4  2  2  9  8  8  8  6  6 10
##  [433]  6  8  6  8 10  5 10  7  8 10 10  6 10  3 10  2  7 10  7  6  7  8  8  7
##  [457]  7  2  6  6 10  7  9  5  7  7  7  1  2 10  8 10  7  9  8  1  8  7  2  8
##  [481]  4  5 10  7 10  6  6  3  9  4  6  5  4 10  2  7  8  7  9  6  2  8 10 10
##  [505]  9 10  7 10  8 10  6  2  7  3  2 10  9 10  8  6  8  9  6  8  2 10  7  4
##  [529]  9  1  8  8  5  1 10  7  8  2  8 10  6 10 10  8  5 10  7  9  2  8 10  6
##  [553]  7 10 10 10  9  2  7 10  9  7  3  9  4  4  6 10  3  9  6  6  6  3  9  7
##  [577]  5  6  9  6  7  1 10  2  5  8  6 10  3  8  4  6  7 10 10 10  4  8  4  8
##  [601]  2  6  8  4  9  8 10  4  8  8  4  5  7 10  7  6  3  7  2  7  8  6  5  8
##  [625]  7  8  7  1  6 10  9  1 10  9  6  7  9  8  9  6  7 10 10 10  2  9  1  8
##  [649] 10  8  9  2  7  8  8 10  9  8  7  9  9  7  8  3  9  7 10 10  2 10  5  9
##  [673]  9 10 10 10  6  7  8 10  1  8  7  9  3  7  7 10  9 10 10  7  8  3  7  2
##  [697] 10  9 10 10  4  4 10  9  6  8 10  8  7  4  3  7  7  4  7  7  8  5  3 10
##  [721]  7 10  6  8  8 10  8  9  6  4  7  3  8  7 10 10  7  4  9  3  1 10  7  2
##  [745]  7  8  8 10 10  5  3  2 10 10  7 10  2 10  7 10  5 10  7  8 10  7  7  6
##  [769]  9 10  8  2 10  5  4  7  7 10 10  4  7  6  8  8  7  7  9  3 10  7  7  9
##  [793]  8  3 10  8  3  8 10  6  2  4 10  7  3 10  9  1  7 10  3 10  7 10 10  2
##  [817] 10  9  7  1  9  2  3 10  7  7 10  7  7  7  5  8 10  7  5 10  7  6  8 10
##  [841]  7  5 10  1  5 10  7 10  7  2 10  7  8 10 10  7  7  8 10  1  1  1 10  9
##  [865]  7  9  2  2  7  9 10 10  7  7  9  9 10  4  7  2  8 10 10 10 10 10  7 10
##  [889]  9  8  3  7  1  8 10 10  7  7  8  1  8  8 10  4  9  5  1  6 10 10  6  6
##  [913]  1  1  9  6  6  7  7  3  4 10  6  8  1  8  8  9 10 10  8  7  7  1  9  8
##  [937] 10 10  5  6  6  7 10 10  5  2  5 10  2  2 10  6 10  4 10  2  2  8  6  8
##  [961] 10  5  8 10  7  8  8  8 10  7 10  3 10  4 10  8  4  6  8 10  6  7  6  3
##  [985]  9  7 10  9  6 10  5  2  7 10  2  9  7  7 10  8

10

1. No passo 2, a distribuição de \(Y\) é Geométrica(\(\lambda\)).

2. O algoritmo simula o processo de gerar uma variável aleatória com uma taxa de hazard discreta constante \(\lambda\). O passo crucial é a verificação em Passo 5. Ele gera múltiplos valores de \(Y\) (que têm distribuição geométrica) e acumula esses valores em \(S\), mas a cada iteração o algoritmo verifica se o processo deve continuar ou terminar com base na distribuição definida pela taxa de hazard. Se a condição for atendida, o valor final de \(Y\) é atribuído a \(X\) e o processo é concluído.

3. Ao gerar \(Y\) com distribuição geométrica e usar a condição de parada no Passo 5 baseada em uma comparação envolvendo \(\lambda\), o algoritmo garante que a probabilidade de sucesso não depende do número de tentativas realizadas até aquele momento. Ou seja, a probabilidade de que \(X\) seja gerado após \(Y\) tentativas segue a mesma lógica de uma distribuição geométrica, com taxa de hazard constante \(\lambda\).

3 Capítulo 5

O capítulo 5 do livro simulation do Sheldon Ross fala sobre a geração de variáveis aleatória contínuas por meio de métodos como a tranformada inversa, aceitação e rejeição, entre outros.

3.1 Questões

1. Dê um método para gerar uma variável aleatória cuja função de densidade é: \[f(x) = \frac{e^x}{e - 1}, \text{ se } 0 \leq x \leq 1\]

2. Dê um método para gerar uma variável aleatória cuja função de densidade é: \[ f(x) = \begin{cases} \frac{x - 2}{2}, & \text{se } 2 \leq x \leq 3 \\ \frac{2 - x/3}{2}, & \text{se } 3 < x \leq 6 \end{cases} \]

3. Use o método da transformada inversa para gerar uma variável aleatória cuja função de distribuição acumulada é: \[F(x) = \frac{x^2 + x}{2}, \text{ se } 0 \leq x \leq 1\]

4. Apresente um método para gerar uma variável aleatória com a função de distribuição \[F(x) = 1 - exp(-\alpha x^\beta), \text{ se } 0 \leq x \leq \infty\] Uma variável aleatória com essa distribuição é chamada de variável aleatória Weibull.

5. Apresente um método para gerar uma variável aleatória com a seguinte função de densidade: \[ f(x) = \begin{cases} e^{2x}, \text{ se } -\infty \leq x \leq 0 \\ e^{-2x}, \text{ se } 0 \leq x \leq \infty \end{cases}\]

6. Apresente dois algoritmos para gerar uma variável aleatória com a seguinte função de distribuição: \[F(x) = 1 - e^{-x} - e^{-2x} + e^{-3x}, \text{ se } x>0\]

7. Apresente dois algoritmos para gerar uma variável aleatória com a seguinte função de densidade: \[f(x) = \frac{1}{4} + 2x^3 + \frac{5}{4}x^4, \text{ se } 0 < x < 1\]

8. Apresente um algoritmo para gerar uma variável aleatória com a seguinte função de densidade: \[f(x) = 2xe^{-x^2}, x>0\]

9. Use o método da rejeição para encontrar uma maneira eficiente de gerar uma variável aleatória com a seguinte função de densidade:\[f(x) = \frac{1}{2}(1 + x)e^{-x}, \text{ se } 0 < x < \infty\]

10. Apresente um método eficiente para gerar uma variável aleatória X com a seguinte função de densidade: \[f(x) = \frac{1}{0,000336}x(1 - x)^3, \text{ se } 0,8 < x < 1\]

3.2 Respostas

1

f1 <- function(x){ # fdp dada na questão
  exp(x)/(exp(1) - 1)
}

g1 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- log(u*(exp(1) - 1) + 1)
}
n <- 100
npa <- g1(n) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X")
curve(f1, add = T, col = "red", lwd = 2)

npa

##   [1] 0.385198316 0.940480140 0.800583211 0.175596308 0.310271406 0.238586606
##   [7] 0.859579036 0.461677837 0.985309546 0.241481714 0.040068944 0.993078191
##  [13] 0.505452047 0.002891558 0.836837295 0.598504769 0.218935564 0.636357490
##  [19] 0.313385749 0.702657596 0.720977452 0.672046416 0.073404844 0.381120589
##  [25] 0.527546108 0.876636963 0.229655233 0.334066641 0.448747536 0.371944859
##  [31] 0.469843304 0.879491755 0.913774263 0.837231375 0.251104420 0.933791987
##  [37] 0.839218966 0.485121716 0.869409265 0.327691240 0.231936371 0.051704965
##  [43] 0.527638820 0.331714141 0.961677363 0.321203858 0.688197059 0.789347490
##  [49] 0.458585500 0.978932588 0.913102532 0.016077722 0.172504774 0.218819092
##  [55] 0.532592760 0.371317956 0.920682264 0.981083169 0.041600857 0.464617689
##  [61] 0.011641210 0.935993760 0.768057404 0.092261664 0.975894809 0.769312272
##  [67] 0.551007810 0.893967912 0.017607747 0.627231439 0.730706806 0.927580927
##  [73] 0.868645190 0.800384592 0.087050151 0.599161981 0.563746711 0.671169524
##  [79] 0.997335654 0.889302995 0.560501742 0.359855689 0.845682947 0.238024893
##  [85] 0.509117310 0.496050445 0.268205883 0.317925182 0.505490744 0.507685341
##  [91] 0.217862619 0.800093219 0.081609597 0.364715803 0.102268079 0.215167430
##  [97] 0.822822600 0.308673012 0.750691606 0.066904472

2

f21 <- function(x){ # primeira parte da fdp dada na questão
  (x-2)/2
}

f22 <- function(x){ # segunda parte da fdp dada na questão
  (2 - (x/3))/2
}

g2 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- numeric(n)
  
  for (i in 1:n){
    if (u[i] <= 1/4){
      x[i] <- 2*sqrt(u[i]) + 2
    }
    else {
      x[i] <- 6 - 2*sqrt(3*(1 - u[i])) 
    }
  }
  return(x)
}
    
n <- 100
npa <- g2(n) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X")
curve(f21, from = 2, to = 3, add = T, col = "red", lwd = 2, ylim = c(0, 0.6), xlim = c(2, 6))
curve(f22, from = 3, to = 6, add = T, col = "red", lwd = 2)

npa

##   [1] 4.614222 3.406421 3.458683 3.331791 4.687352 4.869156 4.700681 4.199780
##   [9] 3.330918 2.114355 3.694665 4.493614 3.726431 3.691756 3.720224 3.395148
##  [17] 4.019495 2.978061 4.773195 3.865537 3.709147 2.648793 3.770790 4.000103
##  [25] 4.990333 3.717016 2.389556 3.381476 3.640166 3.255109 3.019243 3.654884
##  [33] 5.220344 2.986846 5.577280 4.699908 3.608501 4.292067 3.439861 3.140512
##  [41] 4.435620 4.163467 2.966743 5.211741 3.087445 3.913420 5.209402 3.892582
##  [49] 3.014264 2.078590 3.773798 2.844222 3.125850 4.751422 3.026257 4.086052
##  [57] 3.962313 3.466670 4.216172 3.284332 4.846969 3.395931 3.378629 2.652130
##  [65] 2.615916 4.985148 4.450100 5.229668 3.166738 4.798413 3.474706 3.999502
##  [73] 4.340217 3.178604 3.129149 2.644896 2.746053 2.499795 5.318910 4.876938
##  [81] 5.425749 3.868474 2.976908 3.989917 3.179879 4.401246 2.329668 2.614992
##  [89] 4.164570 2.510200 3.148402 3.815494 3.692427 3.313131 4.870582 3.850428
##  [97] 2.671097 3.371180 3.412609 2.961754

3

f3 <- function(x){ # fdp dada na questão
  (2*x + 1)/2
}

g3 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- (sqrt(1 + 8*u) - 1)/2
}
n <- 100
npa <- g3(n); # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X")
curve(f3, add = T, col = "red", lwd = 2)

npa

##   [1] 0.105088002 0.994021751 0.981675346 0.641313685 0.789038154 0.900525828
##   [7] 0.998230214 0.213729269 0.169165289 0.594868907 0.404611075 0.466037511
##  [13] 0.723826933 0.888821962 0.498815787 0.826970160 0.167831750 0.834711716
##  [19] 0.378228566 0.538569954 0.343898199 0.881002402 0.113134996 0.528923439
##  [25] 0.237486969 0.604139181 0.872947939 0.972747556 0.535411460 0.852642235
##  [31] 0.423249450 0.940992701 0.175890672 0.432830950 0.520782660 0.139815934
##  [37] 0.719873864 0.972017994 0.885305875 0.156042344 0.576665903 0.997913012
##  [43] 0.207324908 0.534468581 0.973928494 0.074235456 0.918943775 0.194458951
##  [49] 0.205583818 0.154264570 0.435817800 0.971430463 0.663987574 0.519104341
##  [55] 0.996438523 0.668079902 0.829703006 0.364385290 0.566666668 0.539527597
##  [61] 0.992900274 0.114667897 0.221928354 0.830561412 0.860906736 0.717139576
##  [67] 0.210362649 0.095160715 0.514046928 0.018500402 0.797983462 0.246981785
##  [73] 0.136769102 0.445465319 0.226172628 0.839698576 0.785075959 0.828306092
##  [79] 0.525062575 0.942709138 0.268817946 0.001214142 0.570265824 0.402283657
##  [85] 0.409660768 0.403164398 0.955157736 0.780208183 0.118994748 0.901860355
##  [91] 0.269346506 0.009623160 0.802923673 0.428914953 0.853593873 0.514603552
##  [97] 0.777470811 0.030880748 0.882376340 0.780895977

4

f4 <- function(x, a, b){ # fdp dada na questão
  a*b*(x^(b-1))*exp(-a*(x^b))
}

g4 <- function(n, a, b){ # método da transformada inversa
  u <- runif(n)
  x <- (-(1/a)*log(1 - u))^(1/b)
}
n <- 100
a <- 10
b <- 5
npa <- g4(n, a, b) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X")
curve(f4(x, a, b), add = T, col = "red", lwd = 2)

npa

##   [1] 0.6148701 0.8303975 0.4548639 0.4435391 0.5997637 0.3279569 0.4911659
##   [8] 0.3769831 0.6333398 0.6713092 0.4692308 0.5131925 0.6305407 0.6592147
##  [15] 0.4522374 0.7311355 0.6408340 0.6536097 0.7388307 0.5156375 0.4784782
##  [22] 0.6661914 0.3366795 0.5522506 0.7745184 0.7618112 0.6131397 0.2996715
##  [29] 0.5883196 0.5633998 0.4883801 0.5909035 0.7736996 0.3306774 0.6939429
##  [36] 0.4666587 0.5952189 0.5513206 0.6341110 0.5167045 0.5987992 0.6126731
##  [43] 0.6545976 0.4560236 0.7237953 0.6216451 0.4255699 0.5475460 0.2004476
##  [50] 0.4178770 0.5878103 0.5934978 0.5786403 0.6084987 0.8074545 0.4474122
##  [57] 0.6308867 0.6747714 0.7112265 0.6832070 0.4443734 0.6105295 0.6290045
##  [64] 0.5576105 0.4652808 0.4954902 0.5000025 0.6888867 0.7638267 0.6152785
##  [71] 0.5365962 0.6640240 0.6516204 0.7208120 0.6022737 0.5968385 0.2873863
##  [78] 0.5100547 0.7858225 0.7344019 0.3318457 0.3358841 0.5237325 0.4813996
##  [85] 0.4772544 0.4141876 0.3317259 0.3905530 0.5554268 0.5382627 0.7044881
##  [92] 0.5052246 0.4678306 0.5189209 0.4132897 0.6394400 0.6489532 0.5616085
##  [99] 0.5077156 0.8529108

5

f51 <- function(x){ # primeira parte da fdp dada na questão
  exp(2*x)
}

f52 <- function(x){ # segunda parte da fdp dada na questão
  exp(-2*x)
}

g5 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- numeric(n)
  
  for (i in 1:n){
    if (u[i] < 1/2){
      x[i] <- (1/2)*log(2*u[i])
    }
    else {
      x[i] <- -(1/2)*log(2*(1-u[i])) 
    }
  }
  return(x)
}
    
n <- 100
npa <- g5(n) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X", ylim = c(0, 1))
curve(f51, from = -5, to = 0, add = T, col = "red", lwd = 2)
curve(f52, from = 0, to = 5, add = T, col = "red", lwd = 2)

npa

##   [1] -0.0379797678 -0.0756399701 -0.5894259480 -0.2229259765 -0.8279343839
##   [6] -1.0367352558 -0.3814041928 -1.3689057394 -0.1593211260 -0.6417895844
##  [11] -0.1173596318  0.3861129042 -0.1438074570 -0.2688939642 -0.5646657105
##  [16] -0.9197150190 -2.1971293433 -1.0146906443  0.8481133727 -0.0477841333
##  [21] -1.4240807952 -0.1975869324 -1.9340240513 -1.3151571371  0.1293149545
##  [26]  0.0229214779  0.2467080751  0.5946333187 -1.1041177867 -0.1415330671
##  [31] -0.1483285231  0.2757316762  0.2471003601 -0.2794760146 -0.3617027405
##  [36]  0.3460619377  0.0526998515 -1.2133229482 -0.0568932999  0.2574795272
##  [41]  0.2268012544  0.2617698890 -0.6147695760  0.0450444425  0.9413529321
##  [46]  0.8013116451  0.0190338109  1.0635968170 -0.6425445714 -0.6245793314
##  [51] -0.3843806393  0.1013350321  1.5438940265  0.2554342307 -0.4390661075
##  [56]  0.9332776357 -0.0889792364  0.3770360148 -0.1242874008 -0.7136061583
##  [61] -0.9194245310 -0.0570655919 -0.5492115339  0.0070453258  0.3983534131
##  [66] -0.1924130073  0.0115998827  0.3242181378 -0.6051656881 -0.4030401596
##  [71] -0.1352159220  0.2183671152 -1.2984012698  0.0558122993 -0.4572200824
##  [76]  0.4680156785  0.1691922301  0.4634548307  1.4771188887 -0.2806382422
##  [81] -0.1729241408 -0.8205901122 -1.1821119990  0.0007323967  0.1327448669
##  [86] -0.1917051707 -0.0123141390 -0.6224329814  2.2481611413  0.1770922602
##  [91]  0.5884851769 -0.6025665485  0.2589063160 -0.3586765104  0.5515691874
##  [96] -0.6141788484 -0.0052496404  0.4355199824  0.0902977394  0.0566102898

6

Transformada Inversa

f6_1 <- function(x){ # fdp dada na questão
  exp(-x) + 2*exp(-2*x) - 3*exp(-3*x)
}

g6_1 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- -log(1 - u) # inversa da exponencial com lambda = 1
  
  y <- x + (1/2)*x - (1/3)*x
  return(y)
}
n <- 100
npa <- g6_1(n) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X", ylim = c(0, 0.8))
curve(f6_1(x), add = T, col = "red", lwd = 2)

npa

##   [1] 1.81598691 3.38449057 0.03964653 0.94285451 1.76559030 1.61427824
##   [7] 0.03974705 0.65639901 0.66557002 0.19188301 1.84012919 1.04650122
##  [13] 0.27636548 0.09316581 0.54742124 1.39199648 0.52714103 1.28957854
##  [19] 0.31460720 0.52291718 0.63970977 1.78443329 1.23769520 2.66014081
##  [25] 0.68808627 0.11410693 2.76321595 3.93594639 1.69544482 0.69312189
##  [31] 0.94144784 1.78527175 2.73009865 0.57127022 1.90540038 0.18357927
##  [37] 1.10701143 0.69009439 0.79292833 2.57501124 1.31030270 0.11968624
##  [43] 2.26610956 1.52322414 0.14966136 0.27120881 0.71863160 0.83456859
##  [49] 0.31513305 1.02670561 0.80121951 2.16962005 0.37916629 1.23730714
##  [55] 0.28243549 0.19118401 0.03023575 1.71741830 3.03073287 0.13549618
##  [61] 1.74726862 0.29813872 0.05243160 0.48714149 0.63726391 0.37739507
##  [67] 2.99716837 1.39500346 0.12578256 0.17679580 1.14685786 0.97807781
##  [73] 1.08293637 0.24504727 2.45607450 1.12565357 1.68461995 2.08529990
##  [79] 0.86956238 0.63297935 0.35527143 5.70587805 1.54111311 0.11912451
##  [85] 0.73964576 0.43548113 0.33857086 0.81811192 0.68998957 0.03702394
##  [91] 0.20625636 1.71385473 0.00579918 0.23281157 0.73934522 1.66774905
##  [97] 1.60341055 0.04066956 0.39189150 0.46358688

Aceitação e Rejeição

f6_2 <- function(x){ # fdp dada na questão
  exp(-x) + 2*exp(-2*x) - 3*exp(-3*x)
}

g6_2 <- function(x){ # fdp fácil de gerar
  exp(-x)
}

h6 <- function(n){
  a <- seq(0, 10, by = 0.01)
  c <- max(f6_2(a)/g6_2(a))
  
  i <- 1 #Contador de valores aceitos.
  
  x<-numeric(n)
  
  while(i <= n){ # método da aceitação e rejeição
    u1 <- runif(1)
    y <- -log(1 - u1) # inversa da exponencial 1
    
    u <- runif(1)
    alpha <- f6_2(y)/(c*g6_2(y))
    
    if (u < alpha){
      x[i] <- y
      i <- i + 1
    }
  }
  return(x)
}

n <- 100
npa <- h6(n)

hist(npa, probability = T, main = "Histograma de X", ylim = c(0, 0.8))
curve(f6_2(x), add = T, col = "red", lwd = 2)

npa

##   [1] 0.26004332 0.99205593 1.93179411 0.60426998 0.42467195 1.14233605
##   [7] 0.60788759 0.23331917 0.85922464 1.45619279 0.32020878 0.77562865
##  [13] 2.32500984 1.79354097 2.28251570 0.27405453 1.31177628 0.22895309
##  [19] 1.48516765 0.45230000 1.14225312 2.25040383 0.45867417 0.81847989
##  [25] 0.48802328 1.37940627 3.62408612 1.74089745 0.41275019 1.24461391
##  [31] 1.38446719 0.92801189 2.81029268 0.32356129 2.96264188 0.94719355
##  [37] 1.21129223 0.97020407 1.21658488 2.43968721 1.21304857 1.48600585
##  [43] 1.07856208 0.78279231 0.53951820 0.75726844 1.88898564 0.91202515
##  [49] 1.12297760 0.30188843 0.49307404 0.14702215 0.95368590 0.37610440
##  [55] 2.40690885 1.01243815 0.43251179 0.62536428 1.03260558 0.66486136
##  [61] 0.45409780 0.43788774 0.04619158 0.50852337 1.53287998 0.28805765
##  [67] 0.49615109 1.11766080 1.32547385 0.82783407 0.57591506 1.02064325
##  [73] 1.35361752 0.49448912 0.98695368 0.75665878 0.14728419 1.59839288
##  [79] 1.63266479 0.87147613 0.08655394 0.35915130 1.84533584 0.40180539
##  [85] 0.35219307 1.88931352 1.50259950 0.62442359 1.92371995 1.20850552
##  [91] 2.12216977 0.51300485 0.77676230 0.08024037 4.30192140 0.83847855
##  [97] 1.31715296 0.80699610 0.38057004 2.71084817

7

f7 <- function(x){ # fdp dada na questão
  (1/4) + 2*(x^3) + (5/4)*(x^4)
}

g7 <- function(x){ # fdp fácil de gerar
  1 + 0*x
}

h7 <- function(n){
  a <- seq(0, 1, by = 0.01)
  c <- max(f7(a)/g7(a))
  
  i <- 1 #Contador de valores aceitos.
  
  x <- numeric(n)
  
  while(i <= n){ # método da aceitação e rejeição
    u1 <- runif(1)
    y <- runif(1)
    
    u <- runif(1)
    alpha <- f7(y)/(c*g7(y))
    
    if (u <= alpha){
      x[i] <- y
      i <- i + 1
    }
  }
  return(x)
}

n <- 100
npa <- h7(n)

hist(npa, probability = T, main = "Histograma de X")
curve(f7(x), add = T, col = "red", lwd = 2)

npa

##   [1] 0.99819975 0.94766677 0.55916013 0.99387994 0.67088774 0.89163181
##   [7] 0.60665662 0.55898473 0.54646348 0.99621099 0.98868456 0.93675914
##  [13] 0.55802020 0.57031396 0.86333399 0.60389227 0.64803368 0.97394963
##  [19] 0.43288067 0.82449455 0.90442154 0.17366478 0.76652578 0.91067629
##  [25] 0.08853469 0.57884711 0.94835432 0.23345656 0.82192760 0.90072091
##  [31] 0.67711159 0.65467512 0.79596651 0.92449422 0.71567827 0.83111910
##  [37] 0.30121193 0.90060605 0.96113091 0.81305396 0.84259878 0.75405695
##  [43] 0.60471360 0.92579853 0.28777888 0.61753720 0.86550432 0.97529154
##  [49] 0.91916521 0.44095128 0.62503671 0.86753710 0.41632966 0.72693131
##  [55] 0.94794183 0.84500748 0.55348554 0.90388886 0.66591355 0.91545658
##  [61] 0.96700137 0.79949216 0.95878916 0.89022740 0.97137123 0.52138271
##  [67] 0.98889449 0.42404520 0.65216945 0.99428160 0.73827741 0.91781850
##  [73] 0.52382880 0.94744347 0.91424482 0.98667507 0.53276488 0.77436813
##  [79] 0.75264575 0.95965705 0.89046567 0.79714433 0.05928212 0.51768837
##  [85] 0.70916933 0.54301053 0.71470185 0.84232381 0.79556242 0.81970084
##  [91] 0.82909879 0.71778820 0.88880583 0.88486235 0.71162435 0.81432517
##  [97] 0.90316500 0.95950650 0.73544388 0.55009863

8

f8 <- function(x, a, b){ # fdp dada na questão
  2*x*exp(-(x^2))
}

g8 <- function(n){ # método da transformada inversa
  u <- runif(n)
  x <- (-log(1 - u))^(1/2)
}
n <- 100
npa <- g8(n) # números pseudo aleatórios da fdp dada

hist(npa, probability = T, main = "Histograma de X")
curve(f8(x), add = T, col = "red", lwd = 2)

npa

##   [1] 1.2399000 1.3658524 0.9779983 1.3200842 0.6736670 0.7801067 1.2533167
##   [8] 0.3173513 0.5786644 0.5462868 1.0019120 1.3791939 0.9946356 0.9582260
##  [15] 0.4494747 0.9132237 0.4402564 0.6009189 2.5217422 1.0274674 2.1933433
##  [22] 0.3623285 0.6097498 0.7735349 0.2463886 1.0065155 0.3271649 0.2238604
##  [29] 0.2235241 1.2407765 0.4269852 0.3192296 0.3142015 0.8117155 1.0233221
##  [36] 1.2954633 0.9495955 1.5368138 0.3638602 1.0221101 0.3508170 1.0941104
##  [43] 0.7166457 0.6368697 0.7815770 1.0504408 0.2875667 0.5086939 1.1987613
##  [50] 0.5983639 0.9247495 1.0533128 1.0838104 0.3524509 0.7477692 0.8386263
##  [57] 0.8806353 0.1472645 1.0204599 0.5335743 0.3332418 0.7434480 2.5381600
##  [64] 1.1624600 1.6334412 1.9439312 2.2523983 1.5048289 1.5436992 1.4282930
##  [71] 1.7785061 0.3066705 0.2856022 0.6133894 0.3881637 0.7201889 0.4055001
##  [78] 0.3153075 1.6686902 0.4663282 0.3288854 0.5882480 0.5125396 0.7105960
##  [85] 1.0752151 2.1212162 0.6485075 0.3667761 1.1084448 1.1987417 1.4722595
##  [92] 0.6593554 1.1926778 0.6296544 0.2967160 1.1856572 1.1315371 0.6957084
##  [99] 1.0187867 2.1619333

9

f9 <- function(x){ # fdp dada na questão
  (1/2)*(1+x)*exp(-x)
}

g9 <- function(x){ # fdp fácil de gerar
  exp(-x)
}

h9 <- function(n){
  a <- seq(0, 10, by = 0.01)
  c <- max(f9(a)/g9(a))
  
  i <- 1 #Contador de valores aceitos.
  
  x <- numeric(n)
  
  while(i <= n){ # método da aceitação e rejeição
    u1 <- runif(1)
    y <- -log(1 - u1) # inversa da exponencial(1)
    
    u <- runif(1)
    alpha <- f9(y)/(c*g9(y))
    
    if (u < alpha){
      x[i] <- y
      i <- i + 1
    }
  }
  return(x)
}

n <- 100
npa <- h9(n)

hist(npa, probability = T, main = "Histograma de X")
curve(f9(x), add = T, col = "red", lwd = 2)

npa

##   [1] 1.020708755 3.676436824 0.758451528 2.618115104 1.053309064 2.258665470
##   [7] 1.972387567 3.188518905 3.373589449 0.155932236 0.648409502 0.133054740
##  [13] 0.648100298 1.249811526 1.742375356 2.634290170 1.150555030 0.254271334
##  [19] 0.992829552 0.921346982 0.282017904 0.866198533 0.731046919 3.155588966
##  [25] 1.133462738 0.727170437 1.536468604 0.934014041 1.933589929 1.597662214
##  [31] 6.126614233 2.288763167 2.925913131 1.359458701 4.429854981 4.255398341
##  [37] 3.628831172 2.670500723 0.253149106 0.350618572 0.151696586 1.519148825
##  [43] 0.333295476 0.811710455 0.022723454 1.217749778 0.384869720 3.262775966
##  [49] 0.762261414 1.655633786 0.348738423 1.375660957 1.065756841 1.617702992
##  [55] 2.723250792 1.059110030 0.237904381 1.300225299 0.293507138 2.483514606
##  [61] 1.251246171 0.771365210 1.503671253 0.934117416 0.215857831 0.813612418
##  [67] 0.268146875 2.327592772 0.001314052 1.163566910 2.261681110 2.045382362
##  [73] 1.718973998 0.296574890 2.413745424 0.109499463 0.880697415 0.200747545
##  [79] 0.126963019 0.635045203 0.665649100 0.142527456 2.602725976 3.234055960
##  [85] 0.691999255 0.342740620 1.894095170 2.482407262 3.894867921 1.056495492
##  [91] 1.431741381 0.405597774 3.570399173 2.076314378 2.430593992 2.449925898
##  [97] 3.498794531 3.859100445 6.226008263 0.543308773

10

f10 <- function(x){ # fdp dada na questão
  (1/0.000336)*x*(1-x)^3
}

g10 <- function(x){ # fdp fácil de gerar
  (1/0.2) + 0*x
}

h10 <- function(n){
  a <- seq(0.8, 1, by = 0.01)
  c <- max(f10(a)/g10(a))
  
  i <- 1 #Contador de valores aceitos.
  
  x <- numeric(n)
  
  while(i <= n){ # método da aceitação e rejeição
    u1 <- runif(1)
    y <-  0.2*u1 + 0.8# inversa da uniforme(0.8, 1)
    
    u <- runif(1)
    alpha <- f10(y)/(c*g10(y))
    
    if (u < alpha){
      x[i] <- y
      i <- i + 1
    }
  }
  return(x)
}

n <- 100
npa <- h10(n)

hist(npa, probability = T, main = "Histograma de X")
curve(f10(x), add = T, col = "red", lwd = 2)

npa

##   [1] 0.8997857 0.8375573 0.8356843 0.8317148 0.8441961 0.8236120 0.8496425
##   [8] 0.8700064 0.8466931 0.8374123 0.8032448 0.8896856 0.8095650 0.8831248
##  [15] 0.8968715 0.8129082 0.8314236 0.8089567 0.8686625 0.8473454 0.8714966
##  [22] 0.8281933 0.8274843 0.8732078 0.8110994 0.8219720 0.8856399 0.9353719
##  [29] 0.8706920 0.8247802 0.8026484 0.8029787 0.8099834 0.8062403 0.8577282
##  [36] 0.8587927 0.8308728 0.8644877 0.8942151 0.8166775 0.8910087 0.8336385
##  [43] 0.8011839 0.8620009 0.8067206 0.8220164 0.8328212 0.8124863 0.8137804
##  [50] 0.8126484 0.8795223 0.8048006 0.8358279 0.8036522 0.8866636 0.8887811
##  [57] 0.9008698 0.8449648 0.8160569 0.8119904 0.8644511 0.8120121 0.8188106
##  [64] 0.8346389 0.8128814 0.8410428 0.8163314 0.8083283 0.9038885 0.8922660
##  [71] 0.9504896 0.8269022 0.9318217 0.8241409 0.8152153 0.8477237 0.8570451
##  [78] 0.8250477 0.8276579 0.8417366 0.8817163 0.9392715 0.8294736 0.8350935
##  [85] 0.8334809 0.8385378 0.8203439 0.8584392 0.8176885 0.8543695 0.8338699
##  [92] 0.8050873 0.8280461 0.8497491 0.8118756 0.8276105 0.8095316 0.8143208
##  [99] 0.8487955 0.8748665