Distribuzioni Discrete

Distribuzione Binomiale

Lanciamo una moneta 5 volte. Qual è la probabilità di ottenere 3 teste?

La densità o p(x=x_i) si ottiene con dbinom(x, n, p)

n <- 5          # cinque tentativi
x <- 3          # vogliamo 3 "successi"
p <- 0.5        # una moneta normale con 50% probabilità di testa
dbinom(x, n, p)
## [1] 0.3125

Lanciamo una moneta 5 volte. Qual è la probabilità di ottenere 0,1,..5 teste?

x <- 0:5         # vogliamo 0,1,2,3,4 or 5 "successi"
myProb <- dbinom(x, n, p)
myProb
## [1] 0.03125 0.15625 0.31250 0.31250 0.15625 0.03125
barplot(myProb,names.arg=x)

Discussione

  • Perché la distribuzione è simmetrica? Perché la probabilità di testa o croce è uguale, infatti è completamente arbitrario che chiamiamo testa “successo”.
  • Prima di eseguire in commando: quale sarà il risultato? Perché? Sarà 1, perché la somma delle probabilità deve essere 1.
sum(myProb)
## [1] 1
  • Prima di eseguire in commando: quale sarà il risultato? Perché? Il risultato sarà 0, perché è impossibile di avere 10 successi con 5 tentativi.
dbinom(10, n, p)
## [1] 0


#### Giulio ha scoperto che c’ è una probabilità del 5% che una ragazza accetta un invito a cena. Se chiede a 10 ragazze, qual è la probabilità che nessuna accetta?

n <- 10          # dieci tentativi
x <- 0           # vogliamo 0 "successi"
p <- 0.05        # probabilità di successo 5%
dbinom(x, n, p)
## [1] 0.5987369

Come dipende la probabilità di insuccesso dal numero di tentativi?

Calcolate e producete un grafico.

n <- 0:10        # zero a 10 tentativi
myProb <- dbinom(x, n, p)
plot(n,myProb, main="Probabilita di insuccesso")

Discussione

  • Se aggiungiamo l’argomento type="l" al commando plot, cambia il plot. (Non dimenticate il comma!)
  • Diventa più o meno realistico? Diventa meno realistico perché suggerisce una distribuzione continua.
plot(n,myProb, main="Probabilita di insuccesso",  type="l")

Cambiate il plot così che vada fino a 30 tentativi

n <- 0:30        # zero a 30 tentativi
myProb <- dbinom(x, n, p)
plot(n,myProb,main="Probabilita di insuccesso")

Se prova 10 volte, qual è la probabilità che al massimo 2 ragazze accettano l’ invito?

Versione lenta:

n <- 10          # dieci tentativi
x <- 0:2         # vogliamo 0 a 2 "successi"
dbinom(x, n, p)  # Le probabilità di 0, 1, 2 successi
## [1] 0.5987369 0.3151247 0.0746348
sum(dbinom(x, n, p))  # Prendiamo la somma di queste probabilità
## [1] 0.9884964

Versione veloce:

n <- 10          # dieci tentativi
x <- 2           # probabilità x<=2
pbinom(x, n, p)
## [1] 0.9884964

Discussione

  • Prima di eseguire in commando: quale sarà il risultato? Perché? La probailità di 10 o meno successi con 10 tentativi è 1.
n <- 10               # dieci tentativi
pbinom(10, n, p)
## [1] 1

Se prova 20 volte, qual è la probabilità che almeno 2 ragazze accettano?

n <- 20          # venti tentativi
x <- 1           # vogliamo il complemento della probabilità x<=1
1-pbinom(x, n, p)
## [1] 0.2641605

Alternativa non consigliata: L’ uso della opzione lower.tail.

Se prova 20 volte, può essere al 95% sicuro che il numero di successi sarà inferiore a …

n <- 20          # venti tentativi
q <- 0.95        # quantile del 95%
qbinom(q,n,p)
## [1] 3

Verifica: per x=2, il risultato è sotto il 95% e …

x <- 2
pbinom(x, n, p)
## [1] 0.9245163

… per x=3, il risultato è sopra il 95%.

x <- 3
pbinom(x, n, p)
## [1] 0.9840985

Finalmente vogliamo fare una simulazione.

10 tentativi, 20% probabilità di successo, simuliamo 20 percorsi

n <- 10          # dieci tentativi
p <- 0.2         # probabilità di successo 20%
nSim <- 20       # venti percorsi
rbinom(nSim,n,p)
##  [1] 1 3 2 1 2 2 2 2 2 3 3 0 1 1 3 2 1 2 2 2

Ritorniamo alla moneta del esercizio 1.

Vogliamo simulare 1000 percorsi e produrre un istogramma

n <- 5           # cinque monete
p <- 0.5         # probabilità di successo 50%
nSim <- 1000     # 1000 percorsi
mySim <- rbinom(nSim,n,p)
hist(mySim)      

Non è molto bello … dobbiamo specificare gli intervalli

hist(mySim,breaks=-1:5)

Oppure, alternativamente (più bello)

barplot(table(mySim))

Discussione

  • A che figura di oggi rassomiglia l’ ultima figura? Perché? Rassomiglia alla prima figure in riga 17-21, perché al limite, l’ istogramma converge verso la probabilita

Distribuzione Bernoulli

Non ci esiste una funzione “Bernoulli” in R. Dobbiamo usare quello che abbiamo imparato sulla distribuzione Binomiale: La Bernoulli corrisponde a una Binomiale con n=1.

Qual è la probabilità di avere “testa” con 1 moneta?

n <- 1          # 1 tentativo
x <- 1          # vogliamo 1 "successo"
p <- 0.5        # una moneta normale con 50% probabilità di testa
dbinom(x, n, p)
## [1] 0.5

Discussione

  • Siete sorpresi dal risultato? No!

Adesso simuliamo 100 lanci di una moneta

nSim <- 100     # 100 percorsi
mySim <- rbinom(nSim,n,p)

Qual è il percentuale di “testa” che avete ricevuto?

Ci sono molti possibili soluzioni, per esempio …

mean(mySim)
## [1] 0.39
table(mySim)
## mySim
##  0  1 
## 61 39
summary(mySim)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00    0.39    1.00    1.00
sum(mySim == 1)
## [1] 39

Poisson

In media, passano 50 macchine al minuto sulla A2 Lugano-Mendrisio. Qual è la probabilità che passino esattamente 40 macchine in un certo minuto?

x <- 40         # Siamo interessati in 40 macchine
lambda <- 50    # Il parametro della nostra distribuzione
dpois(x,lambda)
## [1] 0.02149963

Discussione

  • Allora proviamo esattamente 50 macchine in un certo minuto …
  • Prima di eseguire in commando: quale sarà il risultato? Sarà più grande? Perché? Sarà più grande, perché la probabilità è la più grande per \(x=\lambda\).
x <- 50         # Siamo interessati in 50 macchine
dpois(x,lambda)
## [1] 0.05632501

Producete un plot della distribuzione dei numeri di macchine

x <- 0:100
myProb <- dpois(x,lambda)
plot(x, myProb)

Qual è la probabilità che il numero di macchine è al massimo 40?

x <- 40
ppois(x,lambda)
## [1] 0.08607

Qual è la probabilità che il numero di macchine è più grande di 50?

x <- 50          # Più grande di 50 è il complemento di "al massimo 50"
1 - ppois(x,lambda)
## [1] 0.4624833

Qual è la probabilità che passano al massimo 0 macchine? E qual è la probabilità che passano al massimo 100 macchine?

  • Prima di eseguire in commando: quale sarà il risultato? Perché?
ppois(0,lambda)
## [1] 1.92875e-22

Interpretazione: La probabilità è piccola, ma è possibile che non passi nessuna macchina.

ppois(100,lambda)
## [1] 1

Interpretazione: Il risultato “1” suggerisce che che è assolutamente impossibile che passino più che 100 macchine. Questo non e vero. Se proviamo

1 - ppois(100,lambda)
## [1] 1.569745e-10

vediamo che il risultato non è zero. La “1” era un risultato di un arrotondamento

Distribuzioni Continue

Distribuzione uniforme

Se non specifichiamo limiti, R assume la distribuzione uniforme standard tra 0 e 1. La densità della distribuzione uniforme:

dunif(0.5)
## [1] 1
dunif(0.9)
## [1] 1
x <- seq(0,1,0.1)
dunif(x)
##  [1] 1 1 1 1 1 1 1 1 1 1 1

Prima di eseguire in commando: quale sarà il risultato? Perché?

La densita della distribuzione uniforme standard è zero fuori l’ intervallo [0,1].

dunif(1.5)
## [1] 0

Una figura

x <- seq(-0.5,1.5,0.01)
y <- dunif(x)
plot(x,y,type="l")

Per la distribuzione uniforme tra -2 e 2, la densità è

LB <- -2      # LB = "lower bound"
UB <- 2       # UB = "upper bound"
dunif(0,LB,UB)
## [1] 0.25

Adesso producete una figura della densità della distribuzione uniforme tra -1 e 1

x <- seq(-2.5,2.5,0.01)
LB <- -1      # LB = "lower bound"
UB <- 1       # UB = "upper bound"
y <- dunif(x, LB, UB)
plot(x,y,type="l")

Qual è la probabilità che una variable uniforme tra -3 e 3 sia minore di 1?

LB <- -3      # LB = "lower bound"
UB <-  3      # UB = "upper bound"
x  <- 1
punif(x,LB,UB)
## [1] 0.6666667

L’ arrivo di un autobus è distribuito uniforme tra 08:17 è 08:52. Qual è la probabilità che arrivi prima delle ore 08:30?

LB <- 17      # LB = "lower bound"
UB <- 52      # UB = "upper bound"
x  <- 30
punif(x,LB,UB)
## [1] 0.3714286

Gli stipendi in una ditta sono distribuiti uniforme tra CHF 3’500 e 9’770. Possiamo essere sicuri al 90% di non fare più che …

LB <- 3500      # LB = "lower bound"
UB <- 9770      # UB = "upper bound"
q  <- 0.9       # quantile
qunif(q,LB,UB)
## [1] 9143

… e possiamo essere all’ 80% sicuri di fare almeno …

q  <- 0.8
qunif(1-q,LB,UB)
## [1] 4754

Vogliamo simulare i salari di 100 dipendenti di questa ditta

nSim <- 100    # 100 persone
salari <- runif(nSim,LB,UB)

Producete un istogramma dei salari

hist(salari)

Discussione: Rassomiglia una distribuzione uniforme?

Non tanto, perché (a) 100 dipendenti sono pochi, (b) R scegliechiama numeri “rotondi” per i breaks del istogramma. Per un’ istogramma più bello, dobbiamo specificare i breaks:

hist(runif(nSim,LB,UB),breaks=seq(LB,UB,length.out=6))

Qual è la probabilità che una variable con distribuzione uniforme prende esattamente il valore 0.5?

A questa domanda non si può facilmente rispondere con R, ma la risposta e 0 per ogni distribuzione continua.

Distribuzione normale

In assenza di altri parametri, R assume una distribuzione normale standard con \(\mu=0\) e \(\sigma^2=1\). In R, la distribuzione normale si chiama norm

Qual è la densità della distribuzione normale per x=0

x <- 0
dnorm(x)
## [1] 0.3989423

Producete una figura della densità della distribuzione normale

x <- seq(-5,5,0.01)
y <- dnorm(x)
plot(x,y,type="l")

Qual è la probabilità che una variable con distribuzione normale standard sia inferiore di 0?

x <- 0
pnorm(x)
## [1] 0.5

La quantità di pioggia annuale (in mm) a Lugano è distribuita con una distribuzione normale …

… con \(\mu=153\) e \(\sigma=70\). Qual è la probabilità che in un anno la quantità di pioggia sia inferiore di 200mm? NB: R richiede prima la media mu e poi la deviazione standard sd

mu <- 153
sd <- 70
x <- 200
pnorm(x,mu,sd)
## [1] 0.7490262

Qual è la probabilità che in un anno la quantità di pioggia sia superiore di 100mm?

x <- 100
1-pnorm(x,mu,sd)
## [1] 0.7755179

**Qual è la probabilità che in un anno la quantità di pioggia sia tra 100mm e 200mm?

x1 <- 100
x2 <- 200
pnorm(x2,mu,sd)-pnorm(x1,mu,sd)
## [1] 0.5245441

Possiamo essere al 90% sicuri che non piova più di …

q <- 0.9
qnorm(q,mu,sd)
## [1] 242.7086

Il livello massimo del Lago di Lugano ogni anno è distribuito normale con mu=271 e sigma=2

La città vuole costruire una barriera che sostiene una inondazione centenaria. A che altezza deve essere costruita?

mu <- 271
sd <- 2
q = 0.99    # (1 volta in 100 anni)
qnorm(q,mu,sd)
## [1] 275.6527

Finalmente vogliamo simulare il livello massimo del Lago di Lugano per 100 anni

nSim <- 100    # 100 anni
mySim <- rnorm(nSim, mu, sd)

Qual è il valore massimo dei livelli simulati?

max(mySim)
## [1] 275.1057