Pràctica 4: Probabilitat i variables aleatòries discretes

Problema 1: dau piramidal

Un dau piramidal té cinc cares amb els números 2, 3, 5, 7 i 11. El número 11 és a la base de la piràmide i té una probabilitat de sortir que és el doble de la de les altres cares. Considerem la variable aleatòria obtinguda amb els punts que traiem amb una tirada del dau.

  1. Calculeu-ne i representeu-ne la funció de probabilitat.
pesos <- c(1,1,1,1,2)
(fx <- pesos/sum(pesos))
## [1] 0.1666667 0.1666667 0.1666667 0.1666667 0.3333333
x <- c(2, 3, 5, 7, 11)
plot(x, fx, type="h", col="red", ylim=c(0,.5), lwd=3)
points(x, fx, col="red", lwd=3)

  1. Calculeu-ne i representeu-ne la funció de distribució.
(Fx <- cumsum(fx))
## [1] 0.1666667 0.3333333 0.5000000 0.6666667 1.0000000
plot(c(0,x,12), c(0,Fx,1), type="s", col="red", lwd=3)
points(x, Fx, col="red", lwd=3)

  1. Calculeu-ne l’esperança.
(esp <- sum(x*fx))
## [1] 6.5
  1. Calculeu-ne la variància i la desviació estàndard.
(var <- sum(x^2*fx)-esp^2)
## [1] 12.58333
sqrt(var)
## [1] 3.547299
  1. Comproveu els resultats anteriors fent una simulació.
n <- 1e6
tirada1 <- sample(c(2,3,5,7,11), n, replace=TRUE, prob=c(1,1,1,1,2))
table(tirada1)
## tirada1
##      2      3      5      7     11 
## 167104 166603 167301 166541 332451
table(tirada1)/n
## tirada1
##        2        3        5        7       11 
## 0.167104 0.166603 0.167301 0.166541 0.332451
hist(tirada1, breaks=0:12-.5, freq=FALSE)

mean(tirada1)
## [1] 6.49327
var(tirada1)
## [1] 12.57491
sd(tirada1)
## [1] 3.546111
  1. Tirem el dau dos cops. Quina és la probabilitat d’obtenir més de 7 punts sumant les dues tirades?

Indicació: Tot i que és possible calcular analíticament la probabilitat que es demana, fins i tot a mà, es recomana fer una simulació.

# Simulació (resultat aproximat)
n <- 1e7
tirada1 <- sample(c(2,3,5,7,11), n, replace=TRUE, prob=c(1,1,1,1,2))
tirada2 <- sample(c(2,3,5,7,11), n, replace=TRUE, prob=c(1,1,1,1,2))
total <- tirada1+tirada2
table(total>7)/n # Probabilitat de superar el 7
## 
##     FALSE      TRUE 
## 0.1667895 0.8332105
# Podrem també fer un càlcul analìtic (resultat exacte).
# No cal saber-lo fer amb R, tot i que estaria molt bé.
tirades <- c(2,3,5,7,11)
probs <- c(1,1,1,1,2)/6
names(probs) <- tirades # Per etiquetar files i columnes de les matrius
names(tirades) <- tirades
(sumes <- outer(tirades, tirades, "+")) 
##     2  3  5  7 11
## 2   4  5  7  9 13
## 3   5  6  8 10 14
## 5   7  8 10 12 16
## 7   9 10 12 14 18
## 11 13 14 16 18 22
(probsumes <- outer(probs,probs))
##             2          3          5          7         11
## 2  0.02777778 0.02777778 0.02777778 0.02777778 0.05555556
## 3  0.02777778 0.02777778 0.02777778 0.02777778 0.05555556
## 5  0.02777778 0.02777778 0.02777778 0.02777778 0.05555556
## 7  0.02777778 0.02777778 0.02777778 0.02777778 0.05555556
## 11 0.05555556 0.05555556 0.05555556 0.05555556 0.11111111
sumes>7
##        2     3     5    7   11
## 2  FALSE FALSE FALSE TRUE TRUE
## 3  FALSE FALSE  TRUE TRUE TRUE
## 5  FALSE  TRUE  TRUE TRUE TRUE
## 7   TRUE  TRUE  TRUE TRUE TRUE
## 11  TRUE  TRUE  TRUE TRUE TRUE
probsumes[sumes>7]
##  [1] 0.02777778 0.05555556 0.02777778 0.02777778 0.05555556 0.02777778
##  [7] 0.02777778 0.02777778 0.05555556 0.02777778 0.02777778 0.02777778
## [13] 0.02777778 0.05555556 0.05555556 0.05555556 0.05555556 0.05555556
## [19] 0.11111111
sum(probsumes[sumes>7])
## [1] 0.8333333

Problema 2: 1 o 2 daus

En un joc a cada ronda cada jugador tira el dau, i si treu un sis torna a tirar i se sumen els punts de les dues tirades. Considerem la variable aleatòria dels punts obtinguts.

  1. Trobeu-ne la funció de probabilitat i la funció de distribució.
prob <- c(rep(1/6, times=5), 0, rep(1/36, times=6))
x <-1:12
names(prob) <- x          
prob #funció de probabilitat
##          1          2          3          4          5          6 
## 0.16666667 0.16666667 0.16666667 0.16666667 0.16666667 0.00000000 
##          7          8          9         10         11         12 
## 0.02777778 0.02777778 0.02777778 0.02777778 0.02777778 0.02777778
cumsum(prob) #funció de distribució
##         1         2         3         4         5         6         7 
## 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 0.8333333 0.8611111 
##         8         9        10        11        12 
## 0.8888889 0.9166667 0.9444444 0.9722222 1.0000000
  1. Representeu la funció de probabilitat i la de distribució.
fx <- prob
plot(x, fx, type="h", col="red", ylim=c(0,.5), lwd=3)
points(x, fx, col="red", lwd=3)

Fx <- cumsum(fx)
plot(c(0,x,max(x)+1), c(0,Fx,1), type="s", 
     col="red", lwd=3,
     xlab="x")
points(x, Fx, col="red", lwd=3)

  1. Calculeu-ne l’esperança.
(ex <- sum(x*prob))
## [1] 4.083333
  1. Calculeu-ne la variància.
(ex2 <- sum(x^2*prob))
## [1] 24.69444
ex^2
## [1] 16.67361
ex2-ex^2
## [1] 8.020833
  1. Calculeu la probabilitat d’obtenir entre 4 i 7 punts.
sum(fx[x>=4 & x<=7])
## [1] 0.3611111
  1. Comproveu els resultats anteriors amb una simulació.

Problema 3: accidents

En una fàbrica ocorren de mitjana 4,5 accidents al mes.

  1. Quants accidents de mitjana passaran a l’any?

  2. Quina és la probabilitat que passin 5 accidents en un mes? I la que passin 60 accidents en un any?

  3. Quina és la probabilitat que passin més de 3 accidents en un mes? I més de 36 en un any?

  4. Quina és la variància i la desviació estàndard del nombre d’accidents que passa en un any?

  5. Representeu la funció de probabilitat del nombre d’accidents que passa en un any.

  6. Comproveu els resultats anteriors amb una simulació.

Problema 4: cotxes

En un taller d’automoció un 37% dels cotxes que hi arriben tenen una avaria elèctrica. Avui hi entren 9 clients nous amb un cotxe per reparar.

  1. Quants d’aquests cotxes esperarem que tinguin una avaria elèctrica?

  2. Quna és la probabilitat que exactament 3 d’aquests cotxes tinguin una avaria elèctrica?

  3. Quina és la probabilitat que menys de 5 cotxes tinguin una avaria elèctrica.

  4. Calculeu l’esperança i la variància del nombre de cotxes amb una avaria elèctrica.

  5. Representeu la funció de probabilitat de la variable nombre de cotxes que tenen una avaria elèctrica.

  6. Comproveu els resultats anteriors amb una simulació.

Problema 5: Daus

En un país on els jocs d’atzar estan estrictament regulats, es proposa introduir un nou joc on el jugador tira dos daus i el casino en tira tres. Aleshores es compara la suma dels punts dels dos daus del jugador amb la suma dels punts dels tres daus del casino:

  • Si el jugador ha tret més punts que el casino, cobra cinc vegades el que ha apostat.
  • Si el casino ha tret més punts que el jugador, el casino es queda amb el valor de l’aposta.
  • Si empaten l’import de l’aposta va a l’estat en concepte d’impost sobre el joc.

Calculeu, per cada unitat monetària apostada, quant es destina a premis, quant guanya el casino i quant guanya l’estat.

Indicació: feu una simulació. Per calcular les esperances que es demanen caldrà primer saber la probabilitat que guanyi cadascú o que empatin.

# Agafem una mida de la mostra gran. 
n <- 1e7
# Simulem els dos daus del jugador.
jugador <- sample(1:6, size=n, replace = TRUE)+
  sample(1:6, size=n, replace = TRUE)
# Simulem els tres daus del casino.
casino <- sample(1:6, size=n, replace = TRUE)+
  sample(1:6, size=n, replace = TRUE)+
  sample(1:6, size=n, replace = TRUE)
# Estimem les probabilitats comptant els resultats.
(gj <- sum(jugador>casino)/n) #probabilitat de guanyar el jugador
## [1] 0.1521739
(ge <- sum(jugador==casino)/n) #probabilitata de guanyar l'estat
## [1] 0.0693922
sum(jugador<casino)/n #probabilitat de guanyar el casino
## [1] 0.7784339
# I d'aquí surt el destí dels ingressos:
gj * 5 #premis
## [1] 0.7608695
ge #impostos
## [1] 0.0693922
1-ge-5*gj #ingressos del casino
## [1] 0.1697383

Problema 6: Enfilant agulles

Com que tinc poca vista, quan intento enfilar una agulla un 80% de les vegades no encerto el forat i ho he de tornar a intentar.

  1. Quina és la probabilitat que aconsegueixi enfilar una agulla amb cinc intents o menys?
# Dues maneres:
pnbinom(4, 1, .2) #Binomial negativa
## [1] 0.67232
pbinom(0, 5, .2, lower.tail = FALSE) #Bionomial
## [1] 0.67232
  1. He d’enfilar 4 agulles. Quina és la probabilitat que amb 20 intents no ho hagi aconseguit?
pnbinom(20-4, 4, .2, lower.tail = FALSE)
## [1] 0.4114489
# Noteu que també funcionaria:
pnbinom(20.5-4, 4, .2, lower.tail = FALSE)
## [1] 0.4114489
# O
pbinom(3, 20, .2)
## [1] 0.4114489
  1. Quants intents necessito per tal que la probabilitat d’haver enfilat les quatre agulles sigui del 90%?
qnbinom(.9, 4, .2)
## [1] 28

Problema 7: Bicitaxi

La tarifa d’un bicitaxi és de 20 euros els primers 15 minuts i 10 euros cada 15 minuts addicionals. Un 10% dels clients fa viatges fins a 15 minuts, un 25% de 15 a 30 minuts, un 50% de 30 a 45 i la resta de 45 minuts a 1 hora. Considereu la variable aleatòria de l’import d’un viatge.

  1. Trobeu-en la funció de probabilitat i la funció de distribució.
imports <- c(20, 30, 40, 50)
probs <- c(.1, .25, .50, .15)
names(probs) <- imports
probs
##   20   30   40   50 
## 0.10 0.25 0.50 0.15
cumsum(probs)
##   20   30   40   50 
## 0.10 0.35 0.85 1.00
  1. Trobeu-ne l’esperança i la desviació estàndard.
probs*imports
##   20   30   40   50 
##  2.0  7.5 20.0  7.5
(mitjana <- sum(probs*imports))
## [1] 37
imports^2
## [1]  400  900 1600 2500
probs*imports^2
##  20  30  40  50 
##  40 225 800 375
sum(probs*imports^2)
## [1] 1440
mitjana^2
## [1] 1369
(variancia <- sum(probs*imports^2)-mitjana^2)
## [1] 71
sqrt(variancia)
## [1] 8.42615
  1. Trobeu la probabilitat d’obtenir un import més gran de 25 euros.
imports>25
## [1] FALSE  TRUE  TRUE  TRUE
sum(probs[imports>25])
## [1] 0.9
# Alternativa
1-max(cumsum(probs)[imports <= 25])
## [1] 0.9
  1. Comproveu els apartats anteriors amb una simulació.
# a)
n <- 1e6
mostra <- sample(imports, n, replace=TRUE, prob = probs)
table(mostra)
## mostra
##     20     30     40     50 
##  99802 250333 500430 149435
cumsum(table(mostra)/n)
##       20       30       40       50 
## 0.099802 0.350135 0.850565 1.000000
# b)
mean(mostra)
## [1] 36.99498
var(mostra)
## [1] 70.86753
sd(mostra)
## [1] 8.418285
# c)
table(mostra>25)/n
## 
##    FALSE     TRUE 
## 0.099802 0.900198

Pràctica 5: Variables aleatòries contínues

Problema 1: Formigó

La resistència a compressió del formigó que fabrica una central amb una certa dosificació és de 34,2 MPa amb una desviació estàndard de 2,1 MPa.

  1. Quina és la probabilitat que una proveta doni una resistència inferior a 30 MPa?

  2. Quina és la probabilitat que una proveta doni una resistència entre 30 i 35 MPa?

  3. En la denominació del formigó i als càlculs estructurals es fa servir la resistència característica, que és la que és la que té una probabilitat del 95% de ser superada. Quina és la resistència que serà superada pel 95% de les provetes d’aquest formigó?

  4. Representeu la funció de densitat de probabilitat de la resitència del formigó.

  5. Comproveu els resultats anteriors amb una simulació.

Problema 2: oficina d’atenció al ciutadà

El nombre d’usuaris que van cada dia a una certa oficina d’atenció al ciutadà segueix una distribució normal amb mitjana 260 i desviació estàndard 16.

  1. Actualment, el personal de l’oficina té capacitat per atendre 300 usuaris al dia. Quina és la probabilitat que un dia el nombre d’usuaris sigui superior a 300?
pnorm(300, 260, 16, lower.tail = FALSE)
## [1] 0.006209665
  1. Per tal d’estalviar despeses, l’administració decideix reduir el personal de l’oficina de tal manera que la probabilitat de poder atendre tots els usuaris que hi van en un dia sigui del 90%. Quants usuaris han de poder atendre en un dia?
qnorm(.9, 260, 16)
## [1] 280.5048
  1. Després de fer aquest ajust, quina és la probabilitat que puguin atendre tots els usuaris els cinc dies feiners d’una setmana?
dbinom(5, 5, .9)
## [1] 0.59049
# O el que és el mateix:
dbinom(0, 5, .1)
## [1] 0.59049
  1. I la probabilitat que es quedin sense poder atendre tots els usuaris dos dies o més d’una mateixa setmana?
pbinom(3, 5, .9)
## [1] 0.08146
# O el que és el mateix
pbinom(1, 5, .1, lower.tail=FALSE)
## [1] 0.08146

Problema 3: Pi

La simulació també es pot aplicar per resoldre problemes que, a priori, no sembla que tinguin gaire relació amb l’estadística. Un exemple n’és el següent mètode per estimar el valor del nombre \(\pi\).

  1. Considereu un quadrat i el cercle inscrit en aquest quadrat. Si prenem un punt a l’atzar dins del quadrat (amb una distribució uniforme), quina és la probabilitat que caigui dins del cercle?
# La probabilitat és l'àrea del cercle entre l'àrea del quadrat.
# pi/4
  1. Genereu una mostra de punts dins d’un quadrat. Indicació: Genereu una mostra de cada una de les coordenades del punt.
n <- 1e7
x <- runif(n,-1,1)
y <- runif(n,-1,1)
plot(x[1:10000],y[1:10000], pch=".")

  1. Calculeu la distància de cada punt al centre del quadrat i feu-la servir per comptar quina proporció de punts han caigut dins del cercle.
dins <- x^2+y^2<1
sum(dins)/n
## [1] 0.7854119
plot(x[1:10000],y[1:10000],col=c("red","black")[dins+1], pch=".")

  1. Estimeu \(\pi\) a partir de la freqüència anterior.
sum(dins)/n*4
## [1] 3.141648
#Noteu que tot el que hem fet es podria escriure més compacte
sum(runif(n)^2+runif(n)^2<1)/n*4
## [1] 3.140737

Problema 4: Uniforme

Sigui \(X\) una variable uniforme amb mínim=4 i màxim=10.

  1. Representeu-ne la funció de densitat de probabilitat.
x<- seq(3,11, by=.005)
fx<-dunif(x,4,10)
plot(x, fx, type="l", col="red", lwd=3)

Fx<-punif(x,4,10)
plot(x, Fx, type="l", col="red", lwd=3)

  1. Calculeu:
  • \(P(X \leq 6)\)
  • \(P(5 \leq X \leq 8)\)
  • \(P(0 \leq X \leq 20)\)
punif(6, 4,10)
## [1] 0.3333333
punif(8, 4,10)-punif(5, 4,10)
## [1] 0.5
punif(20, 4,10)-punif(0, 4,10) #per aquest no cal ni calcular
## [1] 1
  1. Trobeu la funció de distribució d’aquesta variable.

Problema 5: Tres proves

Una assignatura s’avalua a partir de tres proves: dues proves parcials (25% de la nota cada una) i una prova final (50% de la nota cada una). Si les notes de les tres proves són independents, les notes de les proves parcials es distribueixen uniformement entre 0 i 10 i les de la prova final entre 1 i 10:

  1. Calculeu la probabilitat d’aprovar cada prova.
punif(5, 0, 10, lower.tail=TRUE) #parcial
## [1] 0.5
punif(5, 1, 10, lower.tail=TRUE) #final
## [1] 0.4444444
  1. Calculeu la probabilitat d’aprovar totes les proves.
punif(5, 0, 10, lower.tail=TRUE)^2*punif(5, 1, 10, lower.tail=TRUE)
## [1] 0.1111111
  1. Estimeu la probabilitat d’aprovar l’assignatura. Indicació: Feu una simulació.
n <- 1e6
notes <- 0.25*runif(n, 0, 10) +
  0.25*runif(n, 0, 10) +
  0.5*runif(n, 1, 10)
hist(notes)

sum(notes>=5)/n
## [1] 0.554724

Problema 6: Vendes supermercat

Les vendes diàries d’un supermercat segueixen una distribució normal, amb mitjana 10.000 euros i desviació estàndard 4.000 euros.

Calculeu les probabilitats dels següents esdeveniments:

  1. Que un dia les vendes siguin més petites que 15.000 euros.
mu <- 10000
sigma <- 4000
x <- 15000
(z<- (x-mu)/sigma)
## [1] 1.25
pnorm(x, mu, sigma)
## [1] 0.8943502
  1. Que un dia les vendes siguin més grans de 12.000 euros.
mu <- 10000
sigma <- 4000
x <- 12000
(z<- (x-mu)/sigma)
## [1] 0.5
pnorm(x, mu, sigma, lower.tail = FALSE)
## [1] 0.3085375
  1. Que un dia siguin més grans de 6.000 euros.
mu <- 10000
sigma <- 4000
x <- 6000
(z<- (x-mu)/sigma)
## [1] -1
pnorm(x, mu, sigma, lower.tail = FALSE)
## [1] 0.8413447
  1. Que es trobin entre 10.000 i 15.000 euros.
mu <- 10000
sigma <- 4000
x <- c(15000, 10000)
(z<- (x-mu)/sigma)
## [1] 1.25 0.00
pnorm(x, mu, sigma)
## [1] 0.8943502 0.5000000
pnorm(x, mu, sigma) %*% c(1, -1)
##           [,1]
## [1,] 0.3943502
  1. Que siguin més petites de 7.000 euros.
mu <- 10000
sigma <- 4000
x <- 7000
(z<- (x-mu)/sigma)
## [1] -0.75
pnorm(x, mu, sigma)
## [1] 0.2266274
  1. Quines vendes se superen el 90% dels dies?
mu <- 10000
sigma <- 4000
p <- .9
qnorm(p, lower.tail = FALSE)
## [1] -1.281552
qnorm(p, mu, sigma, lower.tail=FALSE)
## [1] 4873.794
  1. Quines vendes només s’assoleixen un 1% dels dies?
mu <- 10000
sigma <- 4000
p <- .01
qnorm(p, lower.tail=FALSE)
## [1] 2.326348
qnorm(p, mu, sigma, lower.tail=FALSE)
## [1] 19305.39
  1. Si assumim que les vendes de dies diferents són independents, quina distribució segueixen les vendes de dos dies?
# Normal amb:
(mu2d <- mu+mu )
## [1] 20000
(sigma2d <- sqrt(sigma^2+sigma^2))
## [1] 5656.854
  1. Quina és la probabilitat que les vendes de dos dies superin els 24.000 euros?
(z<- (24000-mu2d)/sigma2d)
## [1] 0.7071068
pnorm(24000, mu2d, sigma2d, lower.tail=FALSE)
## [1] 0.2397501

Problema 7: Galetes supermercat

Un 10% dels clients que entren a un supermercat compren galetes. Un dia visiten el supermercat 300 clients.

  1. Quina és la probabilitat que més de 25 clients comprin galetes?
# Amb normal
mu <- .1*300
sigma <- sqrt(300*.1*(1-.1))
x <- 25+.5
(z<- (x-mu)/sigma)
## [1] -0.8660254
pnorm(x, mu, sigma, lower.tail = FALSE)
## [1] 0.8067619
# Amb binomial
pbinom(25, 300, .1, lower.tail = FALSE)
## [1] 0.8050979
  1. I la probabilitat que entre 20 i 40 clients (inclosos) comprin galetes?
# Amb normal
(z <- (c(19.5, 40.5)-mu)/sigma)
## [1] -2.020726  2.020726
pnorm(c(19.5, 40.5), mu, sigma)
## [1] 0.02165407 0.97834593
pnorm(c(19.5, 40.5), mu, sigma) %*% c(-1,1)
##           [,1]
## [1,] 0.9566919
# Amb binomial
pbinom(c(19.5,40.5), 300, .1)
## [1] 0.01711813 0.97461224
  1. Si cada client que compra galetes en compra un sol paquet, quants paquets ha de tenir el supermercat per tenir una probabilitat del 95% de poder tenir galetes per tots els clients que en volen en un dia?
# Amb normal
(z <- pnorm(.95))
## [1] 0.8289439
qnorm(.95, mu, sigma)
## [1] 38.54691
# Amb binomial
qbinom(.95, 300, .1)
## [1] 39

Problema 8: Iogurts

Un dels productes d’una fàbrica de productes làctics és el iogurt amb melmelada de nabius, que es fabrica fent passar l’envàs per una màquina que hi aboca melmedada de nabius i una altra màquina que hi aboca iogurt. La quantitat de producte que hi aboca cada màquina segueix una distribució normal, la del iogurt amb una mitjana de 100 grams i desviació estàndard 5 grams, i la de la melmelada amb una mitjana de 55 grams i una desviació estàndard de 8 grams.

  1. Quina és la probabilitat que el pes net del producte sigui inferior als 150 grams que diu l’etiqueta?

  2. Ens informen que el control de qualitat rebutja el 0,5% de productes que pesin menys. A partir de quin pes els rebutjarà?

Problema 9: Fiabilitat estructural

D’acord amb el projecte una mènsula de formigó té de dimensions nominals de 1,4 m de llum, 1 m d’ample i 18 cm de cantell, però les dimesions reals es distribueixen d’acord amb una distribució normal de mitjana igual a les dimensions nominals i desviació estàndard 1 cm. La densitat del formigó armat també segueix una distribució normal amb mitjana 2400 kg/m3 i una desviació estàndard de 50 kg/m3.

Hem de dimensionar la resistència de la base de la mènsula per tal que resisteixi el seu pes propi. Quin moment de càlcul hem de fer servir per tal que la probabilitat de ser excedit sigui del 1%? I per que sigui de 10-5?

n <- 2e6
llum <- rnorm(n, 1.4, .015)
ample <- rnorm(n, 1.4, .015)
gruix <- rnorm(n, .18, .015)
densitat <- rnorm(n, 2400, 50)
moment <- densitat*ample*gruix*llum/2
hist(moment, breaks="scott")

summary(moment)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   246.5   398.3   423.1   423.4   448.2   606.7
quantile(moment, c(.99, 1-1e-5))
##      99%  99.999% 
## 510.5791 588.5246