#==========================================================================================
# Aula 9 da disciplina ce223 (05/04/2011)
# Distribuição de probabilidades, geração de números aleatórios e estatística descritiva
#                                                               Professor Walmes M. Zeviani
#                                                                     www.leg.ufpr.br/ce223
#==========================================================================================

#------------------------------------------------------------------------------------------
# distribuições discretas de probabilidade

library(gWidgetsRGtk2)

#------------------------------------------------------------------------------------------
# distribuição binomial

par(mfrow=c(1,2))

size <- 10
limits <- list(prob=c(0,1))
plot.dist <- function(...){
  x <- 0:size
  px <- dbinom(x, prob=svalue(prob), size=size)
  pax <- pbinom(x, prob=svalue(prob), size=size)
  plot(x, px, type="h", ylab="P(X=x)", main="Função de probabilidade")
  points(x, px, pch=19)
  abline(v=size*svalue(prob), col=2, lwd=2)
  legend("topleft", c(paste("p =", format(svalue(prob), dig=3)),
                      paste("n =", size)), bty="n")
  legend("topright", c(paste("E(x) =", format(size*svalue(prob), dig=3)),
                       paste("Var(x) =", format(size*svalue(prob)*(1-svalue(prob)), dig=3))), bty="n")
  n <- length(x)
  plot(x, pax, type="n", ylab="P(X<=x)", main="Função de probabilidade acumulada")
  segments(x[-n], pax[-1], x[-1], pax[-1])
  points(x[-n], pax[-1], pch=19)
  points(x[-1], pax[-1], pch=1)
}
w <- gwindow("Controle os parâmetros")
tbl <- glayout(cont=w)
for(i in 1:length(limits)){
  tbl[i,1] <- paste("Deslizador para", names(limits)[i])
  tbl[i,2, expand=TRUE] <- (assign(names(limits)[i],
             gslider(from=limits[[i]][1],
                     to=limits[[i]][2],
                     by=diff(limits[[i]])/30,
                     value=mean(limits[[i]]),
                     container=tbl, handler=plot.dist)))
}
plot.dist()

#------------------------------------------------------------------------------------------
# distribuição Poisson

par(mfrow=c(1,2))

max <- 23
limits <- list(lambda=c(0,30))
plot.dist <- function(...){
  x <- 0:max
  px <- dpois(x, lambda=svalue(lambda))
  pax <- ppois(x, lambda=svalue(lambda))
  plot(x, px, type="h", ylab="P(X=x)", main="Função de probabilidade")
  points(x, px, pch=19)
  abline(v=svalue(lambda), col=2, lwd=2)
  legend("topleft", c(paste("lambda =", format(svalue(lambda), dig=3))), bty="n")
  legend("topright", c(paste("E(x) =", format(svalue(lambda), dig=3)),
                       paste("Var(x) =", format(svalue(lambda), dig=3))), bty="n")
  n <- length(x)
  plot(x, pax, type="n", ylab="P(X<=x)", main="Função de probabilidade acumulada")
  segments(x[-n], pax[-1], x[-1], pax[-1])
  points(x[-n], pax[-1], pch=19)
  points(x[-1], pax[-1], pch=1)
}
source("cria-caixa.R")
plot.dist()

#------------------------------------------------------------------------------------------
# distribições contínuas de probabilidade

#------------------------------------------------------------------------------------------
# distribuição normal

par(mfrow=c(1,2))

range <- c(-4,4)
limits <- list(mean=c(-4,4), sd=c(0.1,3))
plot.dist <- function(...){
  curve(dnorm(x, mean=svalue(mean), sd=svalue(sd)), range[1], range[2],
        ylab="P(X=x)", main="Função densidade de probabilidade")
  abline(v=svalue(mean), col=2, lwd=2)
  legend("topleft", c(paste("mu =", format(svalue(mean), dig=3)),
                      paste("sigma2 =", format(svalue(sd), dig=3))), bty="n")
  legend("topright", c(paste("E(x) =", format(svalue(mean), dig=3)),
                       paste("Var(x) =", format(svalue(sd)^2, dig=3))), bty="n")
  curve(pnorm(x, mean=svalue(mean), sd=svalue(sd)), range[1], range[2],
        ylab="P(X<=x)", main="Função de probabilidade acumulada")
}
source("cria-caixa.R")
plot.dist()

#------------------------------------------------------------------------------------------
# distribuição gamma

par(mfrow=c(1,2))

range <- c(0,40)
limits <- list(shape=c(0,10), scale=c(0.1,5))
plot.dist <- function(...){
  curve(dgamma(x, shape=svalue(shape), scale=svalue(scale)), range[1], range[2],
        ylab="P(X=x)", main="Função densidade de probabilidade")
  abline(v=svalue(shape)*svalue(scale), col=2, lwd=2)
  legend("topleft", c(paste("shape =", format(svalue(shape), dig=3)),
                      paste("scale =", format(svalue(scale), dig=3))), bty="n")
  legend("topright", c(paste("E(x) =", format(svalue(shape)*svalue(scale), dig=3)),
                       paste("Var(x) =", format(svalue(shape)*svalue(scale)^2, dig=3))), bty="n")
  curve(pgamma(x, shape=svalue(shape), scale=svalue(scale)), range[1], range[2],
        ylab="P(X<=x)", main="Função de probabilidade acumulada")
}
source("cria-caixa.R")
plot.dist()

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

rbeta(n, shape1, shape2)
help(rbeta, help_type="html")

#------------------------------------------------------------------------------------------
# distribuição beta

par(mfrow=c(1,2))

range <- c(0,1)
limits <- list(shape1=c(0.01,3), shape2=c(0.01,3))
plot.dist <- function(...){
  curve(dbeta(x, shape1=svalue(shape1), shape2=svalue(shape2)), range[1], range[2],
        ylab="P(X=x)", main="Função densidade de probabilidade")
  abline(v=svalue(shape1)/(svalue(shape1)+svalue(shape2)), col=2, lwd=2)
  legend("topleft", c(paste("shape1 =", format(svalue(shape1), dig=3)),
                      paste("shape2 =", format(svalue(shape2), dig=3))), bty="n")
  legend("topright", c(paste("E(x) =", format(svalue(shape1)/(svalue(shape1)+svalue(shape2)), dig=3)),
                       paste("Var(x) =",
                             format(svalue(shape1)*svalue(shape2)/((svalue(shape1)+svalue(shape2))^2*
                                                                   (svalue(shape1)+svalue(shape2)+1)),
                                    dig=3))), bty="n")
  curve(pbeta(x, shape1=svalue(shape1), shape2=svalue(shape2)), range[1], range[2],
        ylab="P(X<=x)", main="Função de probabilidade acumulada")
}
source("cria-caixa.R")
plot.dist()

#------------------------------------------------------------------------------------------
# geração de números aleatórios

require(fBasics)

#------------------------------------------------------------------------------------------
# gerar uma amostra aleatória de 500 realizações da distribuição normal, mu=3, sd=2.3

aa <- rnorm(500, mean=3, sd=2.3)
str(aa)
basicStats(aa)

#------------------------------------------------------------------------------------------
# como acessar os elementos? como extrair o ic?

bs <- basicStats(aa)
str(bs)

ic <- bs[c("LCL Mean","UCL Mean"),]
ic

#------------------------------------------------------------------------------------------
# como representar graficamente a distribuição dos dados?

layout(1)
par(family="")
hist(aa, freq=FALSE, col=colors()[45])
rug(aa)
lines(density(aa), col=2)
abline(v=mean(aa), col=7, lwd=2)
abline(v=fivenum(aa), col="gray80", lty=2)
abline(v=ic, col="black", lty=2)
curve(dnorm(x, mean=mean(aa), sd=sd(aa)), col=3, add=TRUE)
par(family="mono", font=2)
legend("topleft", capture.output(basicStats(aa)), bty="n")
box()

#------------------------------------------------------------------------------------------
# outro tipo de gráfico

par(family="")
histPlot(as.timeSeries(aa))

#------------------------------------------------------------------------------------------
# gráfico de distribuição acumulada empírica

aa <- rnorm(30, mean=3, sd=2.3)

plot(ecdf(aa))
curve(pnorm(x, mean=3, sd=2.3), col=2, add=TRUE)

plot(ecdf(aa), verticals=TRUE, do.points=FALSE)
curve(pnorm(x, mean=3, sd=2.3), col=2, add=TRUE)

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