Estatística Computacional II

leg.ufpr.br/~walmes/ensino/EC2

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: RColorBrewer
pteo <- 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)