#==========================================================================================
# Aula 10 da disciplina ce223 (08/04/2011)
# Distribuição de probabilidades, integração númerica, esperança e variância
#                                                               Professor Walmes M. Zeviani
#                                                                     www.leg.ufpr.br/ce223
#==========================================================================================

#------------------------------------------------------------------------------------------
# considere a seguinte função

fx <- function(x){ sapply(x, function(xx) ifelse(xx>=0 && xx<=pi/2, sin(xx), 0)) }
fx(seq(0, 5, l=10))

#------------------------------------------------------------------------------------------
# gráfico da função

curve(fx, -pi, pi)

#------------------------------------------------------------------------------------------
# esta é uma função densidade de probabilidade?
# ** satisfaz a condição de ser não negativa.
# ** sua integral é 1?

integrate(fx, 0, pi/2)

** sim, sua integral é 1, portanto umja f.d.p.

#------------------------------------------------------------------------------------------
# qual a esperança variável aleatória representada por fx?

e.fx <- function(x){ x*fx(x) }
class(e.fx)
args(e.fx)
body(e.fx)

Ex <- integrate(e.fx, 0, pi/2)
str(Ex)
Ex$v
abline(v=1)

#------------------------------------------------------------------------------------------
# qual a variância da variável aleatória representada por fx?

v.fx <- function(x){ (x-1)^2*fx(x) }
integrate(v.fx, 0, pi/2)

e.fx2 <- function(x){ x^2*fx(x) }
Ex2 <- integrate(e.fx2, 0, pi/2)
Ex2$v

Ex2$v-Ex$v^2

#------------------------------------------------------------------------------------------
# qual a P(X<1)?

integrate(fx, 0, 1)

#------------------------------------------------------------------------------------------
# qual a P(0.5<x<1.5)

integrate(fx, 0.5, 1.5)

#------------------------------------------------------------------------------------------
# qual a função de distribuição acumulada de X

Fx <- function(x){ 1-cos(x) }
curve(Fx, 0, pi/2)

#------------------------------------------------------------------------------------------
# obter as probs acima agora com a Fx

Fx(1)-Fx(0) # P(0<X<1)
Fx(1.5)-Fx(0.5) # P(0<X<1)

#------------------------------------------------------------------------------------------
# como gerar números aleatórios dessa distribuição?

iFx <- function(u) acos(1-u)
curve(iFx, 0, 1)

aa <- iFx(runif(100000,0,1))
plot(density(aa))
hist(aa, freq=FALSE)
curve(fx, add=TRUE)

#------------------------------------------------------------------------------------------
# veja como funciona

n <- 10
u <- runif(n)
x <- iFx(u)
curve(Fx, 0, pi/2)
rug(x); rug(u, side=2)
segments(rep(0,n), u, x, u, lty=2, col="gray30")
segments(x, rep(0,n), x, u, lty=2, col="gray50")

#------------------------------------------------------------------------------------------
# extra: pintando as áreas correspondentes aos intervalos de integração 0.5 e 1

li <- 0.23; ls <- 1.3
x. <- seq(0,pi/2,l=100)
fx. <- fx(x.)
plot(fx.~x., type="l")
ax <- c(li, li, x.[x.<ls & x.>li], ls, ls)
ay <- c(0,  fx(li), fx.[x.<ls & x.>li], fx(ls), 0)
polygon(ax, ay, dens=12)
integrate(fx, li, ls)

#------------------------------------------------------------------------------------------
# fazer as funções dsen(), psen(), qsen(), e rsen()

dsen <- function(x) fx(x)
dsen(1)

psen <- function(x) Fx(x)
psen(1)

qsen <- function(p) iFx(p)
qsen(0.4596977)

rsen <- function(n) qsen(runif(n))

#------------------------------------------------------------------------------------------
# gerar uma amostra aleatória e obter média e variância

aa <- rsen(1000)
str(aa)

mean(aa)
var(aa)

#------------------------------------------------------------------------------------------