#### MODELLI PROBABILISTICI DISCRETI #####
######################
#### BERNOULLIANA ####
######################
# coins experiment
x <- c()
for(i in 1:30000) x <- c(x,sum(rbinom(100,1,0.5)))
y <- sort(unique(x))
dx <- dbinom(y,100,0.5)
px <- pbinom(y,100,0.5)
plot(y,dx,type="h") # plot(x,dbinom(x,100,0.5))
plot(y,px,type="s") # plot(sort(x),pbinom(sort(x),100,0.5))
lines(ecdf(x),col="green")
var(x) # commentare....
sd(x) # commentare....
# urn experiment (50 balls, 5=white 45=black)
x <- c()
for(i in 1:100) x <- c(x,sum(rbinom(100,1,0.1)))
plot(sort(x),pbinom(sort(x),100,0.1),type="s")
lines(ecdf(x),col="green")
var(x) # commentare....
sd(x) # commentare....
#####################
#### BINOMIALE #####
#####################
# massa di prob (fx)
x <- 0:100
y <- dbinom(x,100,0.2)
plot(x,y,type="h")
points(x,y,pch=20)
title("Distribuzione Binomiale-n=10,p=0.2")
## ripartizione (Fx)
z <- pbinom(x,100,0.2)
plot(x,z,type="s")
points(x,z,pch=20)
## fx (vario p)
par(mfcol=c(3,3))
p <- seq(0.1,0.9,0.1)
col <- rainbow(length(p))
for (i in 1:length(p)) plot(dbinom(0:10,10,p[i]),main=paste("bin n=20, p=",p[i]),type="h",col=col[i],xlab="",ylab="")
## Fx (vario p)
par(mfrow=c(3,3))
p <- seq(0.1,0.9,0.1)
col <- rainbow(length(p))
for (i in 1:length(p)) plot(pbinom(0:1000,1000,p[i]),main=paste("bin n=20, p=",p[i]),type="s",col=col[i],xlab="",ylab="")
# fx(vario n)
p <- 0.2
n <- seq(20,100,10)
col <- rainbow(9)
plot(c(0,dbinom(0:10,10,0.2)),xlim=c(0,30),type="l")
for(i in 1:length(n)) lines(dbinom(0:n[i],n[i],0.2),col=col[i])
# Fx(vario n)
p <- 0.2
n <- seq(20,100,10)
col <- rainbow(9)
plot(pbinom(0:10,10,0.2),xlim=c(0,30),type="l")
for(i in 1:length(n)) lines(pbinom(0:n[i],n[i],0.2),col=col[i])
# qbinom: numero di successi in n prove
# esempio lanciando una moneta 50 volte qual รจ il numero mediano (quantile=50) di successi? qbinom(0.5,50,0.5)
# fare esercizio 5
#########################
#### IPERGEOMETRICA ####
#########################
# formula + esempio urna....
# ipergeometrica e binomiale a confronto
# aumento elementi dell'urna mantenendo costante frazione di successi(p=0.1)
k<-50
m=500
n<-10
x<-c(0:n)
p<-k/m
plot(x,dhyper(x,k,m-k,n),type="b",pch=20)
points(x, dbinom(x,n,p),col="red", pch=20)
lines(x, dbinom(x,n,p),col="red")
k<-500
m=5000
n<-10
x<-c(0:n)
p<-k/m
plot(x,dhyper(x,k,m-k,n),type="b",pch=20)
points(x, dbinom(x,n,p),col="red", pch=20)
lines(x, dbinom(x,n,p),col="red")
k<-5000
m=50000
n<-10
x<-c(0:n)
p<-k/m
plot(x,dhyper(x,k,m-k,n),type="b",pch=20)
points(x, dbinom(x,n,p),col="red", pch=20)
lines(x, dbinom(x,n,p),col="red")
# funzione di ripartizione
k<-50
m=500
n<-10
x<-c(0:n)
p<-k/m
plot(x,phyper(x,k,m-k,n),type="s",ylab="",xlab="",axes=F)
par(new=T)
plot(x,pbinom(x,n,p),type="s",col="red",ylab="",xlab="")
# simulazione dell'estrazione di x elementi da una popolazione ipergeometrica di parametri k, m e n:
k<-20
m<-100
x<-100
n=10
valori_simulati <-rhyper(x, k, m-k, n)
freq_hyper<-table(valori_simulati)/m
plot(freq_hyper, xlim=c(0,n))
# confronto frequenze cumulate teoriche ed empiriche
k<-20
m<-100
n<-10
x<-100000 # rifare la simulazione settando x=10000. cosa succede?
plot(0:n,phyper(0:n,k,m-k,n),type="s")
valori_simulati <-rhyper(x, k, m-k, n)
lines(ecdf(valori_simulati),col="red")
# quantili qhyper(prob, k, m-k,n) lower.tail = TRUE)
k<-200
m=1000
n<-10
prob<-c(0.1,0.25,0.5,0.75,1)
qhyper(prob, k, m-k,n)
######################
#### GEOMETRICA #####
######################
# fx
p <- 0.3
x <- 0:100
plot(dgeom(x,p),pch=20)
points(dgeom(x,0.1),pch=24)
# Fx
p <- 0.3
x <- 0:100
plot(x,pgeom(x,p),type="s",ylim=c(0.1,1))
lines(x,pgeom(x,0.1),type="s",col="red")
# quantili qgeom(p, prob, lower.tail = TRUE)
p<-0.1
prob<-c(0.1,0.25,0.5,0.75,1)
qgeom(prob,p)
# confronto delle frequenze dei valori simulati con grafico il della legge geometrica...cosa succede all'aumentare di m?
p=0.1
m<-10000 # m=1000000
valori_simulati<-rgeom(m,p)
freq_geom<-table(valori_simulati)/m
plot(freq_geom, xlim=c(0,35))
points(0:35, dgeom(0:35,p), col="red", pch=20)
###################
#### POISSON #####
###################
x<-c(0:50)
lambda <-0.6 #0 < lambda <1: funzione decrescente
plot(x, dpois(x, lambda), pch=20)
lambda<-10 #lambda intero: funzione bimodale
plot(x, dpois(x, lambda), pch=20)
plot(x, ppois (x, lambda),type="s")
lambda <-4.5 #lambda non intero>1: funzione unimodale
plot(x, dpois(x, lambda), pch=20)
plot(x, ppois (x, lambda),type="s")
# simulazione dell'estrazione di m elementi da una popolazione poissoniana di parametro lambda e confronto delle frequenze
# dei valori simulati con grafico il della legge di Poisson... cosa succede all'aumentare di m?
lambda <-1
m<-500 # m = 50000
valori_simulati<-rpois(m, lambda)
freq_pois<-table(valori_simulati)/m
plot(freq_pois)
points(0:length(freq_pois), dpois(0:length(freq_pois), lambda),col="red", pch=20)
##################################################################################################################################
##################################################################################################################################
#### MODELLI PROBABILISTICI CONTINUI ###
###################
#### UNIFORME ####
###################
# dunif(x, min=0, max=1, log = FALSE)
dunif(4,3,7)
dunif(1,3,6)
curve(dunif(x,3,7), from=3,to=7)
curve(dunif(x,3,7), from=-2,to=10)
# punif(q, min=0, max=1, lower.tail = TRUE, log.p = FALSE)
curve(punif(x,2,6),-1,10)
curve(punif(x,0,3),-1,10,col="red",add=TRUE)
# quantili qunif(prob, min=0, max=1, lower.tail = TRUE, log.p = FALSE)
prob<-c(0.1,0.25,0.5,0.75,1)
qunif(prob,0,1)
qunif(prob,0,3)
# runif(n, min=0, max=1)
m<-100000 # provare con m=100000
a<-0 # provare anche con a=5, b=10
b<-1
valori_simulati<- runif(m, a, b)
# confronto delle frequenze dei valori simulati con grafico il della legge uniforme
hist(valori_simulati, freq=FALSE)
curve(dunif(x,a,b), xlim=c(a,b), col="red",add=TRUE)
# meglio confrontare le funzioni cumulative, soprattutto trattandosi di variabili casuali continue
plot(ecdf(valori_simulati))
curve(punif(x,a,b), 0, b, add=TRUE ,col="red")
###################
## ESPONENZIALE ##
###################
#dexp(x,rate=1)
dexp(2,4)
dexp(-2,3)
curve(dexp(x,rate=1),0,5)
rate <- c(0.5,2)
for(i in 1:length(rate)) curve(dexp(x, rate=rate[i]),0,5,add=TRUE,col=i+1,lty=i+1)
curve(pexp(x),0,5)
rate <- c(0.5,2)
for(i in 1:length(rate)) curve(pexp(x, rate=rate[i]),0,5,add=TRUE,col=i+1,lty=i+1)
# simulazione...
# rexp(n,rate=1)
m=10000 #poi provare m=10000
rate=0.5
valori_simulati<- rexp(m,rate)
hist(valori_simulati, freq=FALSE, breaks=30)
xinf=5*1/rate
curve(dexp(x,rate), xlim=c(0, xinf),col="red",add=TRUE)
# confronto delle funzioni cumulative
plot(ecdf(valori_simulati))
curve(pexp(x,rate),0, xinf,add=TRUE ,col="red")