Lança histórico da aula 5 e script.

parent 50e01bdb
......@@ -67,9 +67,13 @@ Prof. Dr. [Walmes Zeviani] (walmes [@] ufpr.br)
- Geração de números aleatórios uniformes pelo método da congruência.
4. 12/08: Método da Transformação Integral da Probabilidade.
- Geração de números aleatórios não uniformes pelo método da
transformação integral da probabilidade.
- *Hands on*: programando geração de números da distribuição uniforme e
exponencial.
transformação integral da probabilidade (MTIP).
- *Hands on*: programando geração de números da distribuição uniforme
e exponencial.
5. 17/08: Prática e Soluções Númericas.
- Aplicando o MTIP a partir da função densidade e da acumulada.
- Gerando números da distribuição trapezoidal com interpolação
linear.
## Atividades e avaliações
......
#-----------------------------------------------------------------------
# Gerar números da v.a. X com função de distribuição
# F(x) = x/(\theta + x), x > 0, \theta > 0.
curve(x/(1 + x), from = 0, to = 20)
rx1 <- function(n, theta) {
u <- runif(n)
x <- -u * theta/(u - 1)
return(x)
}
x <- rx1(n = 1000, theta = 1)
# Como saber que gerou corretamente?
plot(ecdf(x))
curve(x/(1 + x), add = TRUE, from = 0, col = 2)
# Ordenar (quantis observados).
x <- sort(x)
# Vetor com as posições dos quantis observados.
i <- seq_along(x)
# Sequência de probabilidades.
p <- (i - 0.5)/length(x)
# Quantis teóricos (para theta = 1).
q <- -p * 1/(p - 1)
# Hum, ruim de ver (escala dominada por extremos).
plot(x ~ q, asp = 1)
# Fazer um log-log quantil-quantil.
plot(x ~ q, asp = 1, log = "xy")
abline(a = 0, b = 1, col = 2)
#-----------------------------------------------------------------------
# Gerar números da v.a. X com função densidade
# f(x) = sin(x)/2, 0 <= x < \pi.
# Gráfico da f(x).
curve(sin(x)/2, 0, pi)
# Integra 1? Verificação numérica.
integrate(function(x) sin(x)/2, 0, pi)
# Gráfico da F(x).
curve(0.5 * (1 - cos(x)), 0, pi)
# Gerador de números aleatórios.
rx2 <- function(n) {
u <- runif(n)
x <- acos(1 - 2 * u)
return(x)
}
x <- rx2(n = 1000)
plot(ecdf(x))
curve(0.5 * (1 - cos(x)), add = TRUE, col = 2)
# Ordenar.
x <- sort(x)
# Vetor com as posições.
i <- seq_along(x)
# Probabilidades.
p <- (i - 0.5)/length(x)
# Quantis teóricos (theta = 1).
q <- acos(1 - 2 * p)
plot(x ~ q, asp = 1)
abline(a = 0, b = 1, col = 2)
#-----------------------------------------------------------------------
# Gerar números aleatórios de X ~ trapeziodal(a, b, c, d).
fx <- function(x, a, b, c, d) {
k <- 2/(d + c - b - a)
y <- 0 +
k * (x - a)/(b - a) * (x >= a & x < b) +
k * (x >= b & x < c) +
k * (d - x)/(d - c) * (x >= c & x < d)
return(y)
}
curve(fx(x, 0, 1, 2, 3), 0, 3)
# Gerar uma sequência de x e avaliar a densidade (y).
b <- 0.1
x <- seq(0, 3, by = b)
y <- fx(x, 0, 1, 2, 3)
plot(y ~ x, type = "o")
points(x, y, type = "h", col = "blue")
n <- length(x)
plot(y ~ x, type = "o")
rect(xleft = x[-n],
ybottom = 0,
xright = x[-1],
ytop = y[-1], col = "gray50")
# Obter a probabilidade.
Px <- cumsum(y * b)
plot(Px ~ x, type = "o")
plot(x ~ Px, type = "o")
# Interpolação linear.
rx3 <- approxfun(x = Px, y = x)
rx3(0.8)
z <- rx3(runif(1000))
plot(ecdf(z))
lines(x = x, y = Px, col = 2)
#-----------------------------------------------------------------------
trapz <- function(x, a, b, c, d) {
k <- 2/(d + c - a - b)
0 +
k * (x - a)/(b - a) * (x >= a & x < b) +
k * (x >= b & x < c) +
k * (d - x)/(d - c) * (x >= c & x < d)
}
curve(trapz(x, 0, 1, 2, 3), 0, 3)
integrate(function(x) trapz(x, 0, 1, 2, 3), 0, 3)
x <- seq(0, 3, by = 0.1)
fx <- trapz(x, 0, 1, 2, 3)
plot(fx ~ x, type = "l")
plot(cumsum(fx * 0.1) ~ x, type = "o")
# Ordenar.
z <- sort(z)
# Vetor com as posições.
i <- seq_along(z)
# Probabilidades.
p <- (i - 0.5)/length(z)
# Quantis teóricos (theta = 1).
q <- rx3(p)
plot(z ~ q, asp = 1)
abline(a = 0, b = 1, col = 2)
#-----------------------------------------------------------------------
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment