Pontificia Universidad Católica de Chile

Facultad de Matemáticas
Departamento de Estadı́stica

Segundo Semestre 2010

EYP1113 - Probabilidad y Estadı́stica

Laboratorio N◦ 5 (Parte 1)

Esta sección del laboratorio tiene como objetivo el calculo de probabilidades y simulación de vari-
ables aleatorias cuyas distribuciones han sido vistas en clases.

Distribución Poisson

Sea X una variable aleatoria con distribución Poisson(λ), es decir, su función de probabilidad y de
distribución acumulada están dada por:
[x]
λk e−λ X λk e−λ
P (X = k) = , P (X ≤ x) =
k! k!
k=0

En R existen funciones implementadas que pueden ser de utilidad:

dpois(x, lambda): P (X = x).

ppois(q, lambda , lower.tail = TRUE) : P (X ≤ q).

ppois(q, lambda , lower.tail = FALSE) : P (X > q).

qpois(p, lambda, lower.tail = TRUE): entrega valor de q tal que P (X ≤ q) = p.

(a) Plante la semilla 1113 para simular una secuencia de números aleatorios entre 0 y 1 de tamaño
10000.

set.seed(1113)
u = runif(10000)

(b) Simule a partir de (a) números aleatorios provenientes de una distribución Poisson de parámetro
λ = 1, 5, 10 y 100.

x1 = qpois(u, lambda = 1)
x2 = qpois(u, lambda = 5)
x3 = qpois(u, lambda = 10)
x4 = qpois(u, lambda = 100)

EYP1113 - Probabilidad y Estadı́stica 1 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
 (c) Utilice la función “table” para obtener la frecuencia relativas de los valores obtenidos. Haga
un “plot” para graficar la frecuencia relativa y utilice la función “lines” para sobreponer la
función de probabilidad teórica.

par(mfrow = c(2,2), cex = 0.8)

plot(table(x1)/10000, ylab = "Frecuencia Relativa", xlab = "x", main = expression(X*" ~ Poisson"*(lambda==1)))
lines(0:6, dpois(0:6, lambda = 1), col = 2, lwd = 2)
abline(h = 0)

plot(table(x2)/10000, ylab = "Frecuencia Relativa", xlab = "x", main = expression(X*" ~ Poisson"*(lambda==5)))

lines(0:15, dpois(0:15, lambda = 5), col = 2, lwd = 2)
abline(h = 0)

plot(table(x3)/10000, ylab = "Frecuencia Relativa", xlab = "x", main = expression(X*" ~ Poisson"*(lambda==10)))

lines(0:23, dpois(0:23, lambda = 10), col = 2, lwd = 2)
abline(h = 0)

plot(table(x4)/10000, ylab = "Frecuencia Relativa", xlab = "x", main = expression(X*" ~ Poisson"*(lambda==100)))

lines(0:140, dpois(0:140, lambda = 100), col = 2, lwd = 2)
abline(h = 0)

X ~ Poisson(λ = 1) X ~ Poisson(λ = 5)
0.15
0.3
Frecuencia Relativa

Frecuencia Relativa

0.10
0.2

0.05
0.1

0.00
0.0

0 1 2 3 4 5 6 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

x x

X ~ Poisson(λ = 10) X ~ Poisson(λ = 100)

0.04
0.12
Frecuencia Relativa

Frecuencia Relativa

0.03
0.08

0.02
0.04

0.01
0.00

0.00

1 2 3 4 5 6 7 8 9 10 12 14 16 18 20 22 64 68 72 76 80 84 88 92 96 101 106 111 116 121 126 132 137

x x

(d) Si X ∼ Poisson(λ = 10), calcule el valor exacto (utilizando dpois o pois) y aproximado (uti-
lizando las frecuencias relativas) de P (5 < X ≤ 15).

Tenemos que P (5 < X ≤ 15) = P (X ≤ 15) − P (X ≤ 5) = P (X = 6) + · · · + P (X = 15).

ppois(15, lambda = 10) - ppois(5, lambda = 10)

0.8841736

sum(dpois(6:15, lambda = 10))

0.8841736

EYP1113 - Probabilidad y Estadı́stica 2 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
 Aproximadamente

(table(x3)/10000)
x3
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0.0003 0.0016 0.0063 0.0201 0.0389 0.0619 0.0906 0.1139 0.1320 0.1194 0.1137 0.0928 0.0698 0.0537 0.0345
16 17 18 19 20 21 22 23
0.0242 0.0134 0.0070 0.0032 0.0014 0.0009 0.0003 0.0001

(table(x3)/10000)[6:15]
x3
6 7 8 9 10 11 12 13 14 15
0.0619 0.0906 0.1139 0.1320 0.1194 0.1137 0.0928 0.0698 0.0537 0.0345

sum((table(x3)/10000)[6:15])
0.8823

Distribución Gamma

Sea X una variable aleatoria con distribución Gamma(α, λ) cuya su función de densidad esta dada
por:
λα α−1 −λ x
fX (x) = x e , x ≥ 0 y α, λ > 0
Γ(α)
En R existen funciones implementadas que pueden ser de utilidad:

dgamma(x, shape = α, rate = λ): fX (x).

pgamma(q, shape = α, rate = λ, lower.tail = TRUE): P (X ≤ q).

pgamma(q, shape = α, rate = λ, lower.tail = FALSE): P (X > q).

qgamma(p, , shape = α, rate = λ, lower.tail = TRUE): entrega valor de q tal que P (X ≤

q) = p.

(e) Simule a partir de (a) números aleatorios provenientes de una distribución Gamma de pará-
metros λ = 1, 5 y α = 1 , 10.

y1 = qgamma(u, shape = 1, rate = 1)

y2 = qgamma(u, shape = 1, rate = 5)
y3 = qgamma(u, shape = 10, rate = 1)
y4 = qgamma(u, shape = 10, rate = 5)

(f) Utilice la función “hist” para obtener el histograma de frecuencia relativas de los valores
simulados en (e) y con la función “lines” sobreponga la función densidad teórica para cada
caso.

y = seq(0,20, 0.01)
par(mfrow = c(2,2), cex = 0.8)
hist(y1, freq = F, main = expression(X*" ~ Gamma"(alpha==1,lambda==1)), xlab = "x", col = "darkgreen")
lines(y, dgamma(y, shape = 1, rate = 1), col = 2, lwd = 2)
hist(y2, freq = F, main = expression(X*" ~ Gamma"(alpha==1,lambda==5)), xlab = "x", col = "darkgreen")
lines(y, dgamma(y, shape = 1, rate = 5), col = 2, lwd = 2)
hist(y3, freq = F, main = expression(X*" ~ Gamma"(alpha==10,lambda==1)), xlab = "x", col = "darkgreen")
lines(y, dgamma(y, shape = 10, rate = 1), col = 2, lwd = 2)
hist(y4, freq = F, main = expression(X*" ~ Gamma"(alpha==10,lambda==5)), xlab = "x", col = "darkgreen")
lines(y, dgamma(y, shape = 10, rate = 5), col = 2, lwd = 2)

EYP1113 - Probabilidad y Estadı́stica 3 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
 X ~ Gamma(α = 1,, λ = 1) X ~ Gamma(α = 1,, λ = 5)

0.8

4
0.6

3
Density

Density
0.4

2
0.2

1
0.0

0
0 2 4 6 8 0.0 0.5 1.0 1.5

x x

X ~ Gamma(α = 10,, λ = 1) X ~ Gamma(α = 10,, λ = 5)

0.12

0.0 0.1 0.2 0.3 0.4 0.5 0.6

0.08
Density

Density
0.04
0.00

5 10 15 20 25 0 1 2 3 4 5

x x

(g) Si X ∼ Gamma(α = 5, λ = 1), calcule el valor exacto (utilizando pgamma) de P (5 < X ≤ 10).

Tenemos que P (5 < X ≤ 10) = P (X ≤ 10) − P (X ≤ 5).

pgamma(10, shape = 5, rate = 1) - pgamma(5, shape = 5, rate = 1)

0.4112406

EYP1113 - Probabilidad y Estadı́stica 4 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
Pontificia Universidad Católica de Chile
Facultad de Matemáticas
Departamento de Estadı́stica

Segundo Semestre 2010

EYP1113 - Probabilidad y Estadı́stica

Laboratorio N◦ 5 (Parte 2)

Esta sección del laboratorio tiene como objetivo estudiar la distribución exacta y aproximada de
funciones de variables aleatorias independientes e idénticamente distribuidas a partir de simula-
ciones.

Suma
2 ) e Y ∼ Normal(µ , σ 2 ) variables aleatorias independientes. En-
Sean X ∼ Normal(µX , σX Y Y
tonces
2
Z = X + Y ∼ Normal(µX + µY , σX + σY2 )

set.seed(2909)
x = rnorm(1000, mean = 10, sd = 4)
y = rnorm(1000, mean = 50, sd = 2)
z = x + y
hist(z, freq = F, col = "gray", main = expression(Z == X+Y * " ~ Normal"(mu[X]+mu[Y],sigma[X]^2+sigma[Y]^2)))
aux = seq(40, 80, 0.01)
lines(aux, dnorm(aux, mean = 10+50, sd = sqrt(4^2+2^2)), lwd = 2, col = 2)

Sean X e Y variables aleatorias independientes con distribución Exponencial de parámetro

λ. Entonces
Z = X + Y ∼ Gamma(2, λ)

x = rexp(1000, rate = 1)
y = rexp(1000, rate = 1)
z = x + y
hist(z, freq = F, col = "gray", main = expression(Z == X+Y * " ~ Gamma"(2,lambda)), ylim = c(0,0.4))
aux = seq(0, 20, 0.01)
lines(aux, dgamma(aux, shape = 2, rate = 1), lwd = 2, col = 2)

Sean X ∼ Poisson(λ) e Y ∼ Poisson(µ) variables aleatorias independientes. Entonces

Z = X + Y ∼ Poisson(λ + µ)

x = rpois(1000, lambda = 2)
y = rpois(1000, lambda = 10)
z = x + y
plot(table(z)/sum(table(z)), freq = F, col = "blue", main = expression(Z == X+Y * " ~ Poisson"(lambda+mu)),
xlim = c(0,30), ylab = "Probability")
abline(h=0)
aux = 0:30
lines(aux, dpois(aux, lambda = 2+10), lwd = 2, col = 2, type = "l", lty = 2)
lines(aux, dpois(aux, lambda = 2+10), lwd = 2, col = 1, type = "h", lty = 2)

EYP1113 - Probabilidad y Estadı́stica 5 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
Teorema Central de Lı́mite

Sean X1 ,. . . , Xn variables aleatorias independientes e identicamente distribuidas con E(Xi ) = µ y

Var(Xi ) = σ 2 . Entonces
n
X aprox.
Xi ∼ Normal(n µ, n σ 2 ),
i=1
cuando n → ∞.
iid
Sean X1 ,. . . , Xn ∼ Exponencial(λ), entonces
n n
X aprox. n
Xi ∼ Gamma(n, λ) ∼ Normal ,
λ λ2
i=1

n = 100
lambda = 5
y = vector("numeric")
m = 1000
for(i in 1:m){
x = rexp(n = 100, rate = lambda)
y[i] = sum(x)
}

hist(y, freq = F, main = expression("Histograma de "*Y==sum(X[i])*", con "*X[i]*" ~ Exponencial"(lambda)))

z = seq(0, 200, 0.01)
lines(z, dgamma(z, rate = lambda, shape = n), col = 2, lwd = 2, lty = 1)
lines(z, dnorm(z, mean = n/lambda, sd = sqrt(n/lambda^2)), col = 4, lwd = 2, lty = 2)

iid
Sean X1 ,. . . , Xn ∼ Uniforme(a, b), entonces
n
(b − a)2
 
X aprox. (a + b)
Xi ∼ Normal n , n
2 12
i=1

n = 100
a = 5
b = 10
y = vector("numeric")
m = 1000
for(i in 1:m){
x = runif(n = n, min = a, max = b)
y[i] = sum(x)
}

hist(y, freq = F, main = expression("Histograma de "*Y==sum(X[i])*", con "*X[i]*" ~ Uniforme"(a,b)), ylim = c(0,0.03))
z = seq(0, 900, 0.01)
lines(z, dnorm(z, mean = n*(a+b)/2, sd = sqrt(n*(b-a)^2/12)), col = 4, lwd = 2, lty = 2)

iid
Sean X1 ,. . . , Xn ∼ Poisson(λ), entonces
n
X aprox.
Xi ∼ Normal (n λ, n λ)
i=1

n = 100
lambda = 1
y = vector("numeric")
m = 1000

EYP1113 - Probabilidad y Estadı́stica 6 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea
 for(i in 1:m){
x = rpois(n = n, lambda = 1)
y[i] = sum(x)
}

plot(table(y)/sum(table(y)), main = expression("Histograma de "*Y==sum(X[i])*", con "*X[i]*" ~ Poisson"(lambda)))

abline(h = 0)
z = 0:200
lines(z, dpois(z, lambda = n*lambda), col = 2, lwd = 2, lty = 1)
lines(z, dnorm(z, mean = n*lambda, sd = sqrt(n*lambda)), col = 4, lwd = 2, lty = 2)

hist(y, freq = F, main = expression("Histograma de "*Y==sum(X[i])*", con "*X[i]*" ~ Poisson"(lambda)))

z = 0:200
lines(z, dpois(z, lambda = n*lambda), col = 2, lwd = 2, lty = 1)
lines(z, dnorm(z, mean = n*lambda, sd = sqrt(n*lambda)), col = 4, lwd = 2, lty = 2)

Mı́nimo y Máximo

Sean X1 , . . . , Xn variables aleatorias continuas independientes e idénticamente distribuidas con fun-

ción de densidad fX y de distribución FX .

Defina
Y1 = mı́n{X1 , . . . , Xn } y Yn = máx{X1 , . . . , Xn }.
Luego
fY1 (y) = n [1 − FX (y)]n−1 fX (y) y fYn (y) = n [FX (y)]n−1 fX (y)

iid
Sean X1 ,. . . , Xn ∼ Exponencial(λ), entonces

Y1 ∼ Exponencial(n λ)

n = 100
lambda = 5
y = vector("numeric")
m = 1000
for(i in 1:m){
x = rexp(n = 100, rate = lambda)
y[i] = min(x)
}

hist(y, freq = F, main = expression("Histograma de "*Y==min*"{"*X[i]*"}"*", con "*X[i]*" ~ Exponencial"(lambda)))

z = seq(0, .2, 0.0001)
lines(z, dexp(z, rate = n*lambda), col = 2, lwd = 2, lty = 1)

EYP1113 - Probabilidad y Estadı́stica 7 Profesores: Ricardo Aravena

Segundo Semestre 2010 Ricardo Olea