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
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)
sum(myProb)
## [1] 1
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
Calcolate e producete un grafico.
n <- 0:10 # zero a 10 tentativi
myProb <- dbinom(x, n, p)
plot(n,myProb, main="Probabilita di insuccesso")
type="l" al commando plot, cambia il plot. (Non dimenticate il comma!)plot(n,myProb, main="Probabilita di insuccesso", type="l")
n <- 0:30 # zero a 30 tentativi
myProb <- dbinom(x, n, p)
plot(n,myProb,main="Probabilita di insuccesso")
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
n <- 10 # dieci tentativi
pbinom(10, n, p)
## [1] 1
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.
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
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
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))
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.
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
nSim <- 100 # 100 percorsi
mySim <- rbinom(nSim,n,p)
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
x <- 40 # Siamo interessati in 40 macchine
lambda <- 50 # Il parametro della nostra distribuzione
dpois(x,lambda)
## [1] 0.02149963
x <- 50 # Siamo interessati in 50 macchine
dpois(x,lambda)
## [1] 0.05632501
x <- 0:100
myProb <- dpois(x,lambda)
plot(x, myProb)
x <- 40
ppois(x,lambda)
## [1] 0.08607
x <- 50 # Più grande di 50 è il complemento di "al massimo 50"
1 - ppois(x,lambda)
## [1] 0.4624833
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
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
La densita della distribuzione uniforme standard è zero fuori l’ intervallo [0,1].
dunif(1.5)
## [1] 0
x <- seq(-0.5,1.5,0.01)
y <- dunif(x)
plot(x,y,type="l")
LB <- -2 # LB = "lower bound"
UB <- 2 # UB = "upper bound"
dunif(0,LB,UB)
## [1] 0.25
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")
LB <- -3 # LB = "lower bound"
UB <- 3 # UB = "upper bound"
x <- 1
punif(x,LB,UB)
## [1] 0.6666667
LB <- 17 # LB = "lower bound"
UB <- 52 # UB = "upper bound"
x <- 30
punif(x,LB,UB)
## [1] 0.3714286
LB <- 3500 # LB = "lower bound"
UB <- 9770 # UB = "upper bound"
q <- 0.9 # quantile
qunif(q,LB,UB)
## [1] 9143
q <- 0.8
qunif(1-q,LB,UB)
## [1] 4754
nSim <- 100 # 100 persone
salari <- runif(nSim,LB,UB)
hist(salari)
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))
A questa domanda non si può facilmente rispondere con R, ma la risposta e 0 per ogni distribuzione continua.
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
x <- 0
dnorm(x)
## [1] 0.3989423
x <- seq(-5,5,0.01)
y <- dnorm(x)
plot(x,y,type="l")
x <- 0
pnorm(x)
## [1] 0.5
… 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
x <- 100
1-pnorm(x,mu,sd)
## [1] 0.7755179
x1 <- 100
x2 <- 200
pnorm(x2,mu,sd)-pnorm(x1,mu,sd)
## [1] 0.5245441
q <- 0.9
qnorm(q,mu,sd)
## [1] 242.7086
mu=271 e sigma=2La 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
nSim <- 100 # 100 anni
mySim <- rnorm(nSim, mu, sd)
max(mySim)
## [1] 275.1057