|
Estatística Computacional II
|
Geração de números aleatórios por relações entre v.a.
Prof. Walmes M. Zeviani
2017-08-16 e 2017-08-18
1 Distribuição Triangular
#-----------------------------------------------------------------------
# Triagular.
# Função densidade de probabilidade.
fx <- function(x) {
ifelse(x < 0, 1 + x, 1 - x) * (-1 < x) * (x < 1)
}
# Gráfico da função.
curve(fx, -1, 1)
#--------------------------------------------
# Gerando números.
# Usando a soma de duas uniformes, i.e. ~U(-1/2, 1/2).
t <- replicate(1000, sum(runif(2, -0.5, 0.5)))
# Poderia-se fazer em um laço também.
# t <- numeric(1000)
# for(i in 1:length(t)) {
# t[i] <- sum(runif(2))
# }
plot(density(t, from = -1, to = 1))
curve(fx, add = TRUE, from = -1, to = 1, col = 2)
# O melhor é verificar na distribuição acumulada.
Fx <- function(x) {
ifelse(x < 0,
x^2/2 + x + 1/2,
-x^2/2 + x + 1/2)
}
plot(ecdf(t))
curve(Fx, add = TRUE, col = 2)
iFx <- function(u) {
ifelse(u < 0.5,
sqrt(2 * u) - 1,
1 - sqrt(2 - 2 * u))
}
curve(iFx, 0, 1)
#--------------------------------------------
# Também pode-se usar os gráficos qq-plot ou pp-plot.
library(lattice)
library(latticeExtra)## Loading required package: RColorBrewerpteo <- ppoints(n = length(t), a = 1/2)
qteo <- iFx(pteo)
qobs <- sort(t)
pobs <- Fx(x = qobs)
# pp-plot.
xyplot(pobs ~ pteo, col = 1, cex = 0.2) +
layer(panel.abline(a = 0, b = 1, col = 2))
# qq-plot.
xyplot(qobs ~ qteo, col = 1, pch = 1, alpha = 0.2)
2 Binomial pela soma de Bernoulli
#--------------------------------------------
# Gerando a soma.
# Valores os parâmetros da binomial
size <- 6
prob <- 0.5
# Números.
s <- 1000
b <- replicate(s, sum(runif(size) > prob))
# Valores teóricos.
x <- 0:size
px <- dbinom(x, size = size, prob = prob)
# Comparação.
cbind(px, freq = table(b)/s)## px freq
## 0 0.015625 0.013
## 1 0.093750 0.107
## 2 0.234375 0.245
## 3 0.312500 0.296
## 4 0.234375 0.224
## 5 0.093750 0.102
## 6 0.015625 0.013# Visualizando no histograma.
hist(b, breaks = c(x, max(x) + 1) - 0.5, freq = FALSE)
points(px ~ x, type = "h", lwd = 3, col = 4)
# Probabilidades acumuladas.
# pbinom(q = x, size = size, prob = prob)
Px <- cumsum(px)
# Gráficos com as acumuladas.
plot(ecdf(b))
points(Px ~ x, type = "s", col = 2)
3 Teorema do limite central para geração de normais
#--------------------------------------------
# Gerando pela soma de uniformes.
s <- 1000
x <- replicate(s, sum(runif(6, -1, 1)))
plot(ecdf(x))
curve(pnorm(x, mean = 0, sd = sqrt(2)), add = TRUE, col = 2)
qqnorm(x)
qqline(x, lty = 2)
Lembre-se que a variância da Uniforme é \((\max - \min)^2/12\). Então a soma tem variância \(6 * (\max - \min)^2/12 = (\max - \min)^2/2\). Com isso, o desvio-padrão é \((\max - \min)/\sqrt{2}\).
Como parametrizar uma distribuição uniforme para produzir números da normal padrão?
#--------------------------------------------
# Usando a soma de Exponenciais.
s <- 1000
r <- 10
n <- 10
x <- replicate(s, sum(rexp(n, rate = r)))
plot(ecdf(x))
curve(pnorm(x,
mean = n * (1/r),
sd = sqrt(n * (1/r^2))),
add = TRUE, col = 2)
Como parametrizar uma exponencial para gerar números de uma normal com média \(\mu\) qualquer mas desvio-padrão unitário?
#--------------------------------------------
# Gerando a partir de uma Poisson.
s <- 1000
l <- 10
n <- 10
x <- replicate(1000, mean(rpois(n, lambda = l)))
plot(ecdf(x))
curve(pnorm(x, mean = l, sd = sqrt(l/n)), add = TRUE, col = 2)
qqnorm(x)
qqline(x, lty = 2)
Usando a média de \(n\) números da Poisson de parâmetro \(\lambda\), como fazer para gerar números de uma normal padrão?
4 Distribuição Cauchy
#--------------------------------------------
# Cauchy padrão é a razão de duas normais padrões.
s <- 1000
x <- rnorm(s)/rnorm(s)
plot(ecdf(x), xlim = c(-5, 5))
curve(pcauchy(x), add = TRUE, col = 2)
5 Poisson pela relação com a Exponencial
#--------------------------------------------
# Números da exponencial: tempo entre eventos.
s <- 10000
e <- rexp(s, rate = 1)
head(e)## [1] 0.14429912 0.47233690 0.01003866 3.02529825 0.96552790 0.86055510# Tempo de ocorrência de cada evento na linha do tempo.
E <- cumsum(e)
head(E)## [1] 0.1442991 0.6166360 0.6266747 3.6519729 4.6175008 5.4780559# Pontos para demarcar o intervalo com pedaços de comprimento 1.
b <- 0:(max(E) + 1)
range(b)## [1] 0 10079# Linha do tempo e divisão em intervalos.
plot(rep(0, min(s, 50)) ~ E[1:min(s, 50)])
abline(h = 0, col = "red")
abline(v = b, col = "gray")
# Números da Poisson.
x1 <- table(cut(E, breaks = b))
freq <- prop.table(table(x1)/sum(x1))
# Números gerados da Poisson.
length(x1)## [1] 10079# Valores da distribuição teórica.
xx <- min(x1):max(x1)
px <- dpois(xx, lambda = 1)
cbind(px, freq)## px freq
## 0 0.3678794412 0.3727552337
## 1 0.3678794412 0.3635281278
## 2 0.1839397206 0.1837483877
## 3 0.0613132402 0.0627046334
## 4 0.0153283100 0.0136918345
## 5 0.0030656620 0.0031749181
## 6 0.0005109437 0.0003968648# Acumuladas.
plot(ecdf(x1))
lines(cumsum(px) ~ xx, type = "s", col = 2)
Trocando de algoritmo. Ao inves de gerar uma exponencial de tamanho \(S = 1000\) vamos trocar por \(k = 1000\) repetições usando uma exponencial de tamanho \(s = 10\). Note que \(S = k \times s\). No entanto, o grau de aproximação da distribuição empírica é pior. Por que isso acontece?
k <- 1000
s <- 10
x2 <- replicate(k, {
E <- cumsum(rexp(s, rate = 1))
b <- 0:(max(E) + 1)
x <- as.vector(table(cut(E, breaks = b)))
return(x)
})
x2 <- unlist(x2)
xx <- min(x2):max(x2)
Px <- ppois(xx, lambda = 1)
plot(ecdf(x2), col = 4)
lines(ecdf(x1), col = 2)
lines(Px ~ xx, type = "s", col = 1)
legend("right",
legend = c("Teórica",
"Empírica com S = 10000",
"Empírica com k * s = 10000"),
col = c(1, 2, 4),
lty = 1,
bty = "n")
6 Algumas implementações da distribuição Binomial
As implementações originais estão disponíveis em FERREIRA (2013) em linguagem Java.
6.1 Relação com a Geométrica
rm(list = objects())
# Geométrica com x = {1, 2, ...}.
rgeometric <- function(p) {
x <- 1
while (runif(1) < p) {
x <- x + 1
u <- runif(1)
}
return(x)
}
# Verifica o gerador de uniformes.
xx <- replicate(1000, rgeometric(p = 0.5))
plot(ecdf(xx))
lines(x = 1:max(xx),
y = pgeom(0:(max(xx) - 1), prob = 0.5),
type = "s",
col = 2)
Atenção: A geométrica do R é \(\Pr(x) = p(1 - p)^x\), que é o número de fracassos até o primeiro sucesso, \(x = {0, 1, ...}\). Mas a geometric() está usando \(\Pr(x) = p(1 - p)^{x - 1}\), que é o número de tentativas até o primeiro sucesso, \(x = {1, 2, ...}\).
# Valores para os parâmetros da binomial.
size <- 30
prob <- 0.5
r <- replicate(10000, {
i <- 1
G <- rgeometric(p = prob)
# G <- rgeom(1, prob = p) + 1
while (G <= size & i <= (size + 1)) {
i <- i + 1
G <- G + rgeometric(p = prob)
# G <- G + (rgeom(n = 1, prob = prob) + 1)
}
x <- i - 1
return(x)
})
plot(ecdf(r))
xx <- 0:max(r)
lines(pbinom(xx, size = size, prob = prob) ~ xx,
type = "s",
col = 2)
6.2 Relação com a Exponencial
rm(list = objects())
# Valores para os parâmetros da binomial
size <- 30
prob <- 0.5
if (prob < 0.5) {
p <- 1 - prob
} else {
p <- prob
}
a <- -log(1 - p)
r <- replicate(1000, {
i <- 1
z <- rexp(1)/(size - i + 1)
while (z <= a & i <= (size + 1)) {
i <- i + 1
z <- z + rexp(1)/(size - i + 1)
}
x <- i - 1
return(x)
})
if (prob < 0.5) {
r <- size - r
}
plot(ecdf(r))
xx <- 0:max(r)
lines(pbinom(xx, size = size, prob = prob) ~ xx,
type = "s",
col = 2)
6.3 Método da transformação integral de probabilidade
rm(list = objects())
# Valores para os parâmetros da binomial.
size <- 30
prob <- 0.3
qrob <- 1 - prob
s <- prob/qrob
a <- (size + 1) * s # vem do (-x + size + 1)/x * prob/(1 - qrob)
r <- replicate(10000, {
x <- 0
px <- qrob^size # Pr(X = 0) = (1 - p)^n
u <- runif(1)
while (u > px) {
u <- u - px
x <- x + 1
px <- (a/x - s) * px
}
return(x)
})
plot(ecdf(r))
xx <- 0:max(r)
lines(pbinom(xx, size = size, prob = prob) ~ xx,
type = "s",
col = 2)
7 Conteúdo complementar
- Univariate Distribution Relationships: http://www.math.wm.edu/~leemis/chart/UDR/UDR.html.
8 Próxima aula
- Método Box-Muller para geração de números da Normal.
- Sabatina.
Referências
FERREIRA, D. F. Estatística computacional em java. Editora UFLA, 2013.
| Estatística Computacional II | Prof. Walmes Zeviani |