O programa R inclui funcionalidade para operações com distribuições de probabilidades. Para cada distribuição existem 4 operações básicas:

• calcula a densidade de probabilidade f(x) no ponto
• calcula a função de probabilidade acumulada F(x) no ponto
• calcula o quantil correspondente a uma dada probabilidade
• retira uma amostra da distribuição

As funções com distribuição normal são implementadas usando o termo norm. Por default as funções assumem a distribuição normal padrão \(N(\mu = 0, \sigma^2 = 1)\).

• dnorm
• pnorm
• qnorm
• rnorm

Note que devemos passar o desvio padrão como argumento e não a variância

# Distribuição normal padrão
# dnorm(-1)
# pnorm(-1)
# qnorm(0.975)
# rnorm(10)

#---------------------
# Argumentos da função
args(rnorm)

## function (n, mean = 0, sd = 1) 
## NULL

Podemos cacular probabilidades

Seja X uma v.a. com distribuição \(N(100, 100)\).

• P[X < 95]

pnorm(95, 100, 10)

## [1] 0.3085375

• P[90 < X < 110]

pnorm(110, 100, 10) - pnorm(90, 100, 10)

## [1] 0.6826895

• P[X > 95]

1 - pnorm(95, 100, 10)

## [1] 0.6914625

pnorm(95, 100, 10, lower = F)

## [1] 0.6914625

Gráficos da N(0,1)

par(mfrow=c(1,2), mar = c(5,5,1,3))
plot(dnorm, -3, 3) # curve(dnorm, -3, 3)
plot(pnorm, -3, 3)

Gráficos da \(N(\mu =100 ,\sigma^2 =64)\)

x <- seq(70, 130, len = 100) 
fx <- dnorm(x, 100, 8) 
plot(x, fx, type = "l") 

Fx <- pnorm(x, 100, 8) 
plot(x, Fx, type = "l")

plot(dnorm, -3, 3, xlab = "valores de X", 
     ylab = "densidade de probabilidade") 
title("Distribuicão Normal\n X ~ N(100, 64)") 

f1 = function(x) dnorm(x, 100, 8)
f2 = function(x) dnorm(x, 90, 8);  f3 = function(x) dnorm(x, 100, 15)

plot(f1, 60, 140, ylab = "f(x)", lwd = 3) 
plot(f2, 60, 140, add = T, col = 2 , lwd = 3, lty = 2) 
plot(f3, 60, 140, add = T, col = 3 , lwd = 3, lty = 4) 
legend(110, 0.05, c("N(100,64)", "N(90,64)", "N(100,225)"), 
       col = 1:3, lwd=3, lty=c(1,2,4))

Curva da distribuicao Normal sobre o histograma

x = rnorm(100, 0, 10)
hist(x, prob=TRUE)
lines(sort(x), dnorm( sort(x), 0, 10) )

Note que devemos fazer o histograma com densidade neste caso.

• A distribuição Binomial inclui a Bernoulli como caso particular;
• Funções: dbinom, rbinom, pbinom, qbinom.

Seja X uma v.a. com distribuição Binomial com n = 10 e p = 0.35.

Inicialmente vemos quais são os argumentos da função dbinom.

args(dbinom)

## function (x, size, prob, log = FALSE) 
## NULL

• O valor x refere-se ao valor que estou calculando a probabiliadade de ocorrer,
  
• size é o número de experimentos Bernoulli e
• p é a probabilidade de sucesso em cada experimento Bernoulli.

• Gráfico da função de probabilidade

x <- 0:10 
fx <- dbinom(x, 10, 0.35) 
plot(x, fx, type = "h" , ylab= "P(X=x)")

• Gráfico da distribuição de probabilidade acumulada

x <- 0:10 
Fx <- pbinom(x, 10, 0.35) 
plot(x, Fx, type = "s", ylab = expression( P(x<=x) ) )

• Calcule \(P(X = 7)\) e \(P(X < 8)\).

dbinom(7, 10, 0.35)

## [1] 0.02120302

pbinom(7, 10, 0.35)

## [1] 0.9951787

sum(dbinom(0:7, 10, 0.35))

## [1] 0.9951787

• Calcule \(P(X \geq 8)\)

1 - pbinom(7, 10, 0.35)

## [1] 0.004821265

pbinom(7, 10, 0.35, lower = F)

## [1] 0.004821265

• Calcule \(P(3 < X\leq 6)\)

pbinom(6, 10, 0.35) - pbinom(3, 10, 0.35)

## [1] 0.4601487

sum(dbinom(4:6, 10, 0.35))

## [1] 0.4601487

• Uniforme contínua: runif, punif, dunif, qunif
• Poisson: rpois, …
• Outras distribuições: dbeta, dcauchy, dchisq, dexp, df, dgamma, dgeom, dlnorm, dmultinom, dnbinom, dt, dweibull

sample:

A função sample() não é restrita à distribuição uniforme discreta, podendo ser usada para sorteios, com ou sem reposição (argumento replace, default sem reposição), com a possibilidade de associar diferentes probabilidades a cada elemento (argumento prob, default probabilidades iguais para os elementos).

sample(1:10, 15, rep = T)

# argumentos
args(sample)

## function (x, size, replace = FALSE, prob = NULL) 
## NULL

# sorteio de 3 numeros entre os inteiros de 0 a 20
sample(0:20, 3)

## [1] 19 10  2

# sorteio de 5 numeros entre os elementos de um certo vetor
x <- c(23, 34, 12, 22, 17, 28, 18, 19, 20, 13, 18) 
sample(x, 5)

## [1] 28 19 13 18 12

# sorteio de 10 numeros entre os possiveis resultados do 
# lancamento de um dado, com reposicao

sample(1:6, 10, rep = T)

##  [1] 1 3 6 3 3 1 4 6 5 2

# idem ao anterior, porem agora com a probabilidade 
# de cada face proporcional ao valor da face.

sample(1:6, 10, prob = 1:6, rep = T) 

##  [1] 5 3 5 4 6 4 4 4 4 6

# prob pode ser pesos, pois a propria função irá normalizar
# para encontrar as probabilidades

Podemos construir uma função que calcula a densidade de probabilidade. Por exemplo:

fn <- function(x) { 
     fx <- (1/sqrt(2 * pi * 100)) * exp((-1/200) * (x - 100)^2) 
     return(fx) 
}

# Grafico com area destacada
x <- seq(70, 130, l = 200) 
fx <- fn(x) 
plot(x, fx, type = "l") 
ax <- c(70, 70, x[x <= 95], 95, 95) 
ay <- c(0, fn(70), fx[x <= 95], fn(95), 0) 
polygon(ax, ay, dens = 10)

Podemos utilizar a função integrate para calcular probabilidades quando se trata de variáveis contínuas.

# Funcao integrate
integrate(fn, -Inf, 95) # Acumulada no ponto 95

## 0.3085375 with absolute error < 2.1e-06

integrate(fn, 90, 110)  # P[90 < X < 110]

## 0.6826895 with absolute error < 7.6e-15

integrate(fn, 95, +Inf) # P[X > 95]

## 0.6914625 with absolute error < 8.1e-05

• A função integrate utiliza um método numérico e por isso informa o erro associado.
• Para gerar amostras de uma distribuição diferente das disponiveis no R é necessário buscar algum pacote ou implementar métodos de geração que fogem ao scopo da disciplina.

Em algumas situações queremos visualizar a função densidade de probabilidade de uma distribuição bidimensional. Para isto, podemos utilizar as seguintes funções:

• contour
• image
• persp

Obs: Note que estas funções podem ser utilizadas em outros contextos e não apenas para avaliar alguma distribuição de probabilidade.

require(mvtnorm) # Normal multivariada
sigma <- matrix(c(1,0.5,0.5,1), ncol=2)
x1 = x2 = seq(-3.5, 3.5, by=0.1)
a = expand.grid(x=x1, y=x2) 
y = matrix(dmvnorm(x=a, mean=c(0,0),sigma=sigma),nrow=length(x1))
par(mar = c(4, 4, 3, 3))
contour(x1, x2, y, lwd=1.5, labcex = 1.2, cex.axis=1.2, cex.lab=1.8)

## Grafico de imagem
image(x1, x2, y, col  = topo.colors(20))

# adicionando pontos de uma amostra gerada
x <- rmvnorm(n=100, mean=c(0,0), sigma=sigma, method="chol")
points(x, pch=16, col = 2)

persp(x1, x2, y, theta = 210, phi = 20, expand = 0.75,d=1, zlab="",
      xlab="nu1",ylab="nu2", shade = 0.5, ticktype = "detailed",r=2)

## no lattice
require(lattice)
g = data.frame(x1 = x1, x2 = x2, y = as.vector(y))
wireframe(y ~ x1*x2, data = g, drape = TRUE, col.regions = heat.colors(150),
          colorkey = TRUE,  screen = list(z = -60, x = -60))