#========================================================================================== # 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) #------------------------------------------------------------------------------------------