Construir uma função que gera amostra aleatórias (parâmetro de tamanho n) de uma distribuição de probabilidade bi-triangular conforme discutido em sala.
A função é dada por:
f(x) = 8x, se x ∈ [0,1/4); −8x + 4, se x ∈ [ 1/4 , 1/2 ); 8x − 4, se x ∈ [ 1/2, 3/4); −8x + 8, se x ∈ [ 3/4, 1].
Integrando a função para saber se a mesma é uma f.d.p..
# A soma deve ser igual a 1, para que seja uma fdp
integrate(function(x) 8*x,0,1/4)
## 0.25 with absolute error < 2.8e-15
integrate(function(x) -8*x + 4,1/4,1/2)
## 0.25 with absolute error < 2.8e-15
integrate(function(x) 8*x - 4,1/2,3/4)
## 0.25 with absolute error < 2.8e-15
integrate(function(x) -8*x + 8,3/4,1)
## 0.25 with absolute error < 2.8e-15
Conforme comprovado acima temos que a função f(x) é uma f.d.p.. A partir disto escolhemos o método da inversa generalizada a obrenção do propósito desta atividade.Para utilizar esta técnica é necessário a F(x) da distribuição. A F(x) é dada por:
F(x) = 4x^2, se x ∈ [0, 1/4); -4x^2 + 4x −1/ ,se x ∈ [1/4, 1/2); 4x^2 - 4x +3/2, se x ∈ [1/2, 3/4); −4x2 + 8x − 3, se x ∈ [3/4, 1].
# Função F(x)
plot.new()
plot.window(range(0:1), range(0:1))
axis(1, cex=.5)
axis(2, cex=.5)
mtext("F(x)", side=2, line=3, las=0, cex=1)
mtext("x", side=1, line=3, las=0, cex=1)
curve(4*x^2,0,1/4,add = T)
curve(-4*x^2 + 4*x - 0.5,1/4,1/2,add = T)
curve(4*x^2 - 4*x + (3/4) + 0.75,1/2,3/4,add = T)
curve(-4*x^2 + 8*x - 3,3/4,1,add = T)
abline(h = 1,lty = 2)
Para utilizar o algoritmo da inversa generalizada, é necessario a inversa da F(x). Abaixo, uma possível função em R para gerar números aleatórios para esta distribuição é dado abaixo:
# retu = retornar u
generate_bitri <- function(n,retu = FALSE){
u <- runif(n,0,1)
z <- ifelse(u<0.25,sqrt(u)/2,
ifelse(u<0.5,0.5 - 0.5*sqrt(0.5-u),
ifelse(u<0.75,0.5 + 0.5*sqrt(u - 0.5),
1 - 0.5*sqrt(1- u))))
if(retu==FALSE){return(z)}
else{return(data.frame(z,u))}
}
z <- generate_bitri(1000,retu = TRUE)
## Warning in sqrt(0.5 - u): NaNs produced
## Warning in sqrt(u - 0.5): NaNs produced
# Histograma
hist(z$z, freq = FALSE)
curve(8*x,0,1/4,add = T)
curve(-8*x + 4,1/4,1/2,add = T)
curve(8*x - 4,1/2,3/4, add = T)
curve(-8*x + 8,3/4,1,add = T)
# F(x)
plot.new()
plot.window(range(0:1), range(0:1))
axis(1, cex=.5)
axis(2, cex=.5)
mtext("F(x)", side=2, line=3, las=0, cex=1)
mtext("x", side=1, line=3, las=0, cex=1)
curve(4*x^2,0,1/4,add = T)
curve(-4*x^2 + 4*x - 0.5,1/4,1/2,add = T)
curve(4*x^2 - 4*x + (3/4) + 0.75,1/2,3/4,add = T)
curve(-4*x^2 + 8*x - 3,3/4,1,add = T)
abline(h = 1,lty = 2)
points(z$z,z$u, col = "darkred",pch = 20)
Observa-se que o gerador é adequado devido o comportamento da curva no gráfico.