########################################################################## # Laboratorio di "Metodi Statistici avanzati" con R # Prof.ssa Lea Petrella # Facolta' di Economia # Dipartimento Metodi e Modelli per l'Economia il Territorio e la Finanza # Sapienza Universita' di Roma # a.a. 2015-2016 # alberto.diiorio@uniroma1.it ########################################################################## ########################################################################## # Lezione 3: VARIABILI ALEATORIE ########################################################################## # # Indice: # # * distribuzioni e simulazione di variabili aleatorie discrete e continue # * simulazione di variabili aleatorie discrete: il comando sample # * schema generale per le distribuzioni parameteriche celebri # * variabile aleatoria di Bernoulli di parametro "p" # * variabile aleatoria binomiale di parametri "n" e "p" # * variabile aleatoria di Poisson di parametro "lambda" # * variabile aleatoria geometrica di parametro "p" # * variabile aleatoria binomiale negativa di parametri "k" e "p" # * variabile aleatoria ipergeometrica di parametri "n","h" e "r" # * variabile aleatoria normale di parametri "mu" e "sigma" # * variabile aleatoria chi quadrato con "k" gdl # * variabile aleatoria t di Student con "k" gdl # * variabile aleatoria uniforme nell'intervallo [a,b] # * variabile aleatoria beta di parametri "a" e "b" # * variabile aleatoria F di Fisher con "n" e "m" gdl # * variabile aleatoria esponenziale di parametro "mu" # * variabile aleatoria gamma di parametri "alpha" e "beta"# ########################################################################## #----------------------------------------------------------------- # * simulazione di variabili aleatorie discrete: il comando sample #----------------------------------------------------------------- ?sample # ESEMPIO: # Possiamo creare un'urna contenente # quattro palline bianche numerate da 1 a 4 e tre # nere numerate da 1 a 3: urna = c("b1","b2","b3","b4","n1","n2","n3") urna # ed estrarre due palline senza reinserimento sample(urna, size = 2) # ed estrarre due palline con reinserimento sample(urna, 2, replace = T) # ...e se "size" > del numero degli elementi nell'urna? length(urna) sample(urna, size = 10) # OPS!! sample(urna, 10, replace = T) # Nota: ogni volta che simuliamo otteniamo un campione diverso sample(urna, 10, replace = T) # Col comando set.seed PRIMA del comando di simulazione # possiamo fissare il seme di generazione e riottenere lo stesso # campione set.seed(123) sample(urna, 10, replace = T) set.seed(123) sample(urna, 10, replace = T) # ESEMPIO: # il *NUMERO* di teste per due lanci # di una moneta (bilanciata) puo' essere simulato cosi': # 1° lancio C # 2° lancio C k=0 con p=0.5*0.5=0.25 # 1° lancio C # 2° lancio T k=1 con p=0.5*0.5=0.25 # 1° lancio T # 2° lancio C k=1 con p=0.5*0.5=0.25 # 1° lancio T # 2° lancio T k=2 con p=0.5*0.5=0.25 # dove k = {Numero di teste in 2 lanci} # puo' assumere solo tre valori: 0, 1 o 2. # quinidi # Pr(k=0)=Pr(k=2)=0.25 # Pr(k=1)=0.5 # ...facciamolo in R... k = 0:2 p = c(1,2,1)/4 # guardiamo un po' come e' fatta la distribuzione discreta... plot(k,p,type="h",main="Distr. della v.a. numero di teste in due lanci successivi \n di una moneta regolare",ylim=c(0,1),xlab="Numero di teste",ylab="Probabilita' ") sample(k, size = 1, prob = p) sample(k, size = 10, replace=T, prob = p) # Se estraggo un campione casuale di ampiezza elevata # la distribuzione empirica (cioe' le frequenze relative # osservate nel campione) approssima quella teorica # (cioe' le probabilita' Pr(k=0)=Pr(k=2)=0.25 e Pr(k=1)=0.5) n.teste=sample(k, size = 1000, replace=T, prob = p) plot(table(n.teste)/1000,ylim=c(0,1),xlab="Numero di teste",ylab="Frequenza relativa",main="Distr. empirica del numero di teste \n nel campione simulato") #----------------------------------------------------------------- # * schema generale per le distribuzioni parameteriche celebri #----------------------------------------------------------------- #------------------------------------------- # Per le principali variabili aleatorie, # abbiamo a disposizione 4 funzioni con iniziale specifica: # 1) funzione di densita' # [iniziale - "d"] # 2) funzione di ripartizione # [iniziale - "p"] # 3) quantili # [iniziale - "q"] # 4) generatore di numeri pseudo-casuali # [iniziale - "r"] #------------------------------------------- #----------------------------------------------------------------- # * variabile aleatoria di Bernoulli di parametro "p" #----------------------------------------------------------------- # X e' l'evento risultato di una prova con solo due valori # Supporto di X e': {0,1} # E(X) = 0*(1-p) + 1*p = p # Var(X) = (0-p)^2*(1-p) + (1-p)^2*p = p(1-p) # Per la v.a. di Bernoulli non esistono i comandi specifici! # Possiamo usare "sample" per generare # campioni da questa distribuzione n = 10 p = 1/4 sample(0:1,size = n, replace = TRUE, prob = c(1-p,p)) #----------------------------------------------------------------- # * variabile aleatoria binomiale di parametri "n" e "p" #----------------------------------------------------------------- # K e' l'evento numero di successi in "n" prove Bernoulliane i.i.d. # ove la probabilita' di successo in ogni prova e' p # Supporto di K e': {0,1,2,3,...,n} # E(K) = n*p # Var(K) = n*p(1-p) # Prima di tutto fissiamo i parametri n=2 p=0.5 # esempio fatto all'inizio... # Il coefficiente binomiale che compare nella # definizione di questa distribuzione conta # il numero di modi in cui "k" oggetti possono # essere scelti da "n" oggetti distinti. # Per calcolarlo in R possiamo usare la funzione choose help(choose) # Calcoliamo le prob. relative a ciacun valore del supporto # in corrispondenza dei valori fissati dei parametri # P(K=0) prob.di.0=choose(n,0)*p^0*(1-p)^(n-0) prob.di.0 # P(K=1) prob.di.1=choose(n,1)*p^1*(1-p)^(n-1) prob.di.1 # P(K=2) prob.di.2=choose(n,2)*p^2*(1-p)^(n-2) prob.di.2 # piu' facilmente utilizzando la funzione dbinom dbinom(0,size=2,prob=p) dbinom(1,2,p) dbinom(2,2,p) # il comando calcola le probabilita' anche su piu' # valori contemporaneamente dbinom(0:2,2,p) # NOTA: se n=1 ---> v.a. di Bernoulli # ESEMPIO: # supponiamo di aver lanciato una moneta bilanciata 10 volte. # Allora K ~ Binomiale(10,0.5). # La probabilita' che K = 5 possiamo calcolarla a mano: p=0.5 n=10 choose(n,5)*p^5*(1-p)^(n-5) # Piu' semplicemente dbinom(5, size = n, prob = p) # oppure in un colpo solo calcoliamo tutte le prob dbinom(0:10, size = n, prob = p) # arrotondiamo round(dbinom(0:10, size = n, prob = p),3) # Bene, possiamo ora rappresentare graficamente la distr. di probabilità # della v.a. Binomiale(10,0.5) plot(0:10,dbinom(0:10, size = n, prob = p),xlab="K",ylab="Probabilita' ", main="Distribuzione di prob. della v.a. Binomiale \n di parametri n=10 e p=5",type="h",lwd=2) # FUNZIONE DI RIPARTIZIONE: P(K<=k) # utilizzando "pbinom" # La probabilita' che ci siano NON PIU' di 6 successi: P(K <= 6) # ad esempio e' pari a sum(dbinom(0:6, size = n, prob=p)) # oppure direttamente con "pbinom" pbinom(6, size = n, prob=p) # Se stiamo cercando la probabilita' di ALMENO 7 successi: P(K > 6) sum(dbinom(7:n, size = n, prob = p)) 1 - pbinom(6, n, p ) pbinom(6, size = n, prob = p, lower.tail = FALSE) # Rappresentiamo graficamente la funzione di ripartizione al variare di p n=20 curve(pbinom(x,n,0.1),xlim=c(0,20),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. Binomiale \n al variare del parametro p") curve(pbinom(x,n,0.2),col="blue",add=T,lwd=2) curve(pbinom(x,n,0.3),col="red",add=T,lwd=2) curve(pbinom(x,n,0.4),col="green",add=T,lwd=2) curve(pbinom(x,n,0.5),col="violet",add=T,lwd=2) legend(15,0.9,legend=paste("p",seq(.1,.5,.1),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione binomiale di parametro "p" n=10 p=0.5 qq=seq(0.25,0.75,length=3) qbinom(qq,n,p) # cambiamo valore a p p=0.7 qbinom(qq,n,p) # in generale e' possibile ottenere un qualsiasi quantile # come il primo e nono decile (cioe' a livello del 10% e 90%) qbinom(c(.10,.90),n,p) # o anche il quinto e novantacinquesimo percentile (cioe' a livello del 5% e 95%) qbinom(c(0.05,0.95),n,p) # simulare 10 valori da binom(n,p) rbinom(10,n,p) # Maggiore e' il numero dei campioni simulati, migliore e' l'approssimazione! par(mfrow=c(1,3)) plot(0:n,dbinom(0:10,n,p),main="Distr. di prob. di una v.a. Binomiale \n con parametri n=10 e p=0.5",type="h",lwd=2) N1=100 camp=rbinom(N1,n,p) plot(table(camp)/N1,main="Campione di ampiezza 100 da una v.a. Binomiale \n con parametri n=10 e p=0.5", col="red") N2=10000 camp.grande=rbinom(N2,n,p) plot(table(camp.grande)/N2,main="Campione di ampiezza 10.000 da una v.a. Binomiale \n con parametri n=10 e p=0.5",col="blue") # verifichiamo empiricamente che la media sia N*p n=10 p=0.5 # --> E(K)=5 mean(rbinom(1000,n,p)) #Verifichiamo empiricamente l'approssimazione normale della binomiale quando n è "grande" #ossia Binom(n,p)~N(np,np(1-p)) n=100 p=0.7 camp.TLC=rbinom(N2,n,p) plot(table(camp.TLC/N2,main="Distribuzione empirica di una v.a. Binomiale \n con parametri n=100 p=0.7",col="blue") curve(dnorm(x,70,(100*0.7*0.3)^0.5),xlim=c(53,86),add=T,col="red",lwd=3) #----------------------------------------------------------------- # * variabile aleatoria di Poisson di parametro "lambda" #----------------------------------------------------------------- # X numero di eventi che si registrano in certo intervallo # Supporto di X e': {0,1,2,...} # E(X)=lambda # Var(X)=lambda lambda=4 # P(X=0) prob.di.0=exp(-lambda)*lambda^0/factorial(0) prob.di.0 # P(X=1) prob.di.1=exp(-lambda)*lambda^1/factorial(1) prob.di.1 # P(X=2) prob.di.2=exp(-lambda)*lambda^2/factorial(2) prob.di.2 # piu' facilmente utilizzando la funzione dpois dpois(0,lambda=lambda) dpois(1,lambda) dpois(2,lambda) # o tutte insieme dpois(0:2,lambda) # Rappresentazione grafica della distrb di prob. sui valori fino a 20 # (la v.a. di Poisson ha infatti supporto infinito!) plot(0:20,dpois(0:20,lambda),main="Distr. di prob. di una v.a. di Poisson \n con parametro lambda=4", type="h",xlab="X",ylab="Probabilità ") # FUNZIONE DI RIPARTIZIONE # P(X <= 6) # utilizzando "dpois" sum(dpois(0:6, lambda)) # oppure "ppois" ppois(6, lambda) # P(X < 6) ppois(5, lambda) # P(X > 8) ppois(8, lambda,lower.tail=F) # oppure 1-ppois(8, lambda) # Disegniamo la funz. di ripartizione per diverse configurazioni parametriche curve(ppois(x,1),xlim=c(0,15),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. di Poisson \n al variare del parametro lambda") curve(ppois(x,2),col="blue",add=T,lwd=2) curve(ppois(x,3),col="red",add=T,lwd=2) curve(ppois(x,4),col="green",add=T,lwd=2) curve(ppois(x,5),col="violet",add=T,lwd=2) legend(11,0.9,legend=paste("lambda",1:5,sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione di poisson lambda=4 qq=seq(0.25,0.75,length=3) qpois(qq,lambda) # quantile a livello dell'80% (ovvero 0.8) qpois(.8,lambda) # simulare 10 valori da pois(lambda) rpois(10,lambda) # verifichiamo empiricamente che la media sia lambda lambda=0.2 # --> E(K)=0.2 mean(rpois(10,lambda)) mean(rpois(1000,lambda)) mean(rpois(100000,lambda)) #----------------------------------------------------------------- # * variabile aleatoria geometrica di parametro "p" #----------------------------------------------------------------- # T e' l'evento tempo di primo successo # p e' la probabilita' di successo in ogni prova # Supporto di T e': {1,2,3,...} # P(T=k) = p * (1-p)^(k-1) # NOTA: tra i valori inclusi c'e' anche lo zero (riscalamento): T --> V=T-1 # Supporto di V e': {0,1,2,3,...} # P(V = k) = P(T = k+1) = p * (1-p)^(k) # V e' l'evento numero di insuccessi prima del primo successo # E(V)=(1-p)/p # Var(V)=(1-p)/p^2 p=0.5 # P(V=0) prob.di.0=p*(1-p)^(0) prob.di.0 # P(V=1) prob.di.1=p*(1-p)^(1) prob.di.1 # P(V=2) prob.di.2=p*(1-p)^(2) prob.di.2 #...e cosi' via... # piu' facilmente utilizzando la funzione dgeom dgeom(x=0,prob=p) dgeom(1,p) dgeom(2,p) # FUNZIONE di RIPARTIZIONE # P(V <= 3) = P(T <= 4) # utilizzando "dgeom" sum(dgeom(0:3, prob = p)) # utilizzando "pgeom" pgeom(3,p) # analogamente P(V > 2) = P(T >3) pgeom(2, prob = p, lower.tail=F) #oppure 1-pgeom(2,p) curve(pgeom(x,0.1),xlim=c(0,20),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. Geometrica \n al variare del parametro p") curve(pgeom(x,0.2),col="blue",add=T,lwd=2) curve(pgeom(x,0.3),col="red",add=T,lwd=2) curve(pgeom(x,0.4),col="green",add=T,lwd=2) curve(pgeom(x,0.5),col="violet",add=T,lwd=2) legend(16,0.8,legend=paste("p",seq(0.1,0.5,.1),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione geometrica di parametro "p" qq=seq(0.25,0.75,length=3) qgeom(qq,p) # cambiamo valore a p... p=0.1 qgeom(qq,p) # in generale e' possibile ottenere un qualsiasi quantile qgeom(0.54,p) # generare valori "pseudo-causuali" da una v.a. geometrica di parametro "p" a=rgeom(10000,p) a par(mfrow=c(1,2)) plot(0:max(a),dgeom(0:max(a),p),main="Distribuzione geometrica con p=0.1",type="h",xlab="V") plot(table(a)/10000,main="10000 estrazioni da una v.a. geometrica con p=0.1",type="h",xlab="V",col="red") # verifichiamo empiricamente che la media sia (1-p)/p p=0.5 # --> E(V)=1 mean(rgeom(100000,p)) p=0.2 # --> E(V)=0.8/0.2 = 4 mean(rgeom(100000,p)) #----------------------------------------------------------------- # * variabile aleatoria binomiale negativa di parametri "k" e "p" #----------------------------------------------------------------- # N e' l'evento numero di prove bernoulliane i.i.d. # necessarie per ottenere k successi # Supporto di N e': {k,k+1,k+2,...} # p e' la probabilita' di successo in ogni prova # come nel caso della v.a. geometrica si ha un riscalamento X=N-k ---> # X e' l'evento numero di insuccessi precendeti il successo k-esimo # Supporto di X e': {0,1,2,...} # E(X)=k*(1-p)/p # Var(X)=k*(1-p)/p^2 p=0.5 k=5 # P(X=0) prob.di.0=choose((0+k-1),0)*p^k*(1-p)^(0) prob.di.0 # P(X=1) prob.di.1=choose((1+k-1),1)*p^k*(1-p)^(1) prob.di.1 # P(X=2) prob.di.2=choose((2+k-1),2)*p^k*(1-p)^(2) prob.di.2 # piu' facilmente utilizzando la funzione dnbinom dnbinom(0,size=5,prob=p) dnbinom(1,5,p) dnbinom(2,5,p) # FUNZIONE DI RIPARTIZIONE # P(X <= 6) # utilizzando "dnbinom" sum(dnbinom(0:6, size = k, prob = p)) # oppure "pnbinom" pnbinom(6, size = k, prob = p) # P(X > 6) 1 - pnbinom(6, size = k, prob = p) pnbinom(6, size = k, prob = p, lower.tail = FALSE) # se vogliamo approssimare... sum(dnbinom(7:10000, size = k, prob = p)) k=4 curve(pnbinom(x,k,0.1),xlim=c(0,20),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. Binomiale Negativa \n al variare del parametro p") curve(pnbinom(x,k,0.2),col="blue",add=T,lwd=2) curve(pnbinom(x,k,0.3),col="red",add=T,lwd=2) curve(pnbinom(x,k,0.4),col="green",add=T,lwd=2) curve(pnbinom(x,k,0.5),col="violet",add=T,lwd=2) legend(0,1,legend=paste("p",seq(0.1,0.5,0.1),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione binomiale negativa di parametro "p" p=0.5 k=5 qq=seq(0.25,0.75,length=3) qnbinom(qq,n,p) # cambiamo valore a p... p=0.7 qnbinom(qq,k,p) # in generale e' possibile ottenere un qualsiasi quantile qnbinom(0.95,k,p) # simulare 10 valori da negbinom(k,p) rnbinom(10,k,p) # verifichiamo empiricamente che la media sia k*(1-p)/p k=5 p=0.5 # --> E(X)=5 mean(rnbinom(100000,k,p)) p=0.2 # --> E(K)=5*0.8/0.2=20 mean(rnbinom(100000,k,p)) #----------------------------------------------------------------- # * variabile aleatoria ipergeometrica di parametri "N","h" e "r" #----------------------------------------------------------------- # urna contenente N palline: h bianche e N-h nere # X e' l'evento numero di palline bianche estratte in r estrazioni # (SENZA REIMMISSIONE) # Supporto di X e': max {0, r − (N − h)} ≤ X ≤ min {h, r} # E(X)=r*h/N # Var(X)=r*h/N * (N-h)/N * (N-r)/(N-1) # nota R vuole come argomenti h, N-h e r N=20 h=7 r=5 # P(X=0) prob.di.0=choose(h,0)*choose((N-h),(r-0))/choose(N,r) prob.di.0 # P(X=1) prob.di.1=choose(h,1)*choose((N-h),(r-1))/choose(N,r) prob.di.1 # P(X=2) prob.di.2=choose(h,2)*choose((N-h),(r-2))/choose(N,r) prob.di.2 # piu' facilmente utilizzando la funzione dhyper dhyper(0,h,(N-h),r) dhyper(1,h,(N-h),r) dhyper(2,h,(N-h),r) # FUNZIONE DI RIPARTIZIONE # P(X <= 3) # utilizzando "dhyper" sum(dhyper(0:3, h,(N-h),r)) # oppure "phyper" phyper(3, h,(N-h),r) # P(X > 4) 1 - phyper(4, h,(N-h),r) phyper(4, h,(N-h),r, lower.tail = FALSE) sum(dhyper(5:max(h,r), h,(N-h),r)) # approssimazione... # quartili della distribuzione ipergemoetrica di parametri h,(N-h),r qq=seq(0.25,0.75,length=3) qhyper(qq,h,(N-h),r) # simulare 10 valori da v.a. ipergemoetrica di parametri h,(N-h),r rhyper(10,h,(N-h),r) # verifichiamo empiricamente che la media sia E(X)=r*h/N N=20 h=7 r=5 # --> E(X)=5*7/20=1.75 mean(rhyper(100000,h,(N-h),r)) r=10 # --> E(X)=10*7/20=3.5 mean(rhyper(100000,h,(N-h),r)) #----------------------------------------------------------------- # * variabile aleatoria normale di parametri "mu" e "sigma^2" #----------------------------------------------------------------- # X v.a. continua con supporto R # mu in R , sigma>=0 # E(X) = mu # Var(X) = sigma^2 # di default normale standard N(0,1) # nota in R riparametrizzazione con sigma (sd) e non sigma^2 # funzione di densita' nel punto x=0 x=0 1/sqrt(2*pi)*exp(-x^2/2) # oppure attraverso la funzione dnorm() dnorm(x) # diamo altri valori a mu e sigma mu=1 sigma=3 # --> Var(X)=9 # funzione di densita' nel punto x=0 1/(sqrt(2*pi)*sigma)*exp(-((x-mu)^2/(2*sigma^2))) dnorm(0,mean=1,sd=3) # vediamo come varia la densita' normale al variare dei parametri # al variare di mu (sigma costante e pari a 2) curve(dnorm(x,0,2),xlim=c(-10,10),ylim=c(0,0.25),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. normale \n al variare del parametro mu (sigma=2)") curve(dnorm(x,2,2),col="blue",add=T,lwd=2) curve(dnorm(x,4,2),col="red",add=T,lwd=2) curve(dnorm(x,-2,2),col="green",add=T,lwd=2) curve(dnorm(x,-4,2),col="violet",add=T,lwd=2) legend(6,0.25,legend=paste("mu",c(0,2,4,-2,-4),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # al variare di sigma (mu costante e pari a 0) curve(dnorm(x,0,1),xlim=c(-12,12),ylim=c(0,0.4),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. normale \n al variare del parametro sigma (mu=0)") curve(dnorm(x,0,2),col="blue",add=T,lwd=2) curve(dnorm(x,0,3),col="red",add=T,lwd=2) curve(dnorm(x,0,4),col="green",add=T,lwd=2) curve(dnorm(x,0,5),col="violet",add=T,lwd=2) legend(5,0.4,legend=paste("sigma",1:5,sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # FUNZIONE DI RIPARTIZIONE # Data la seguente configurazione di parametri mu=2 sigma=3 # calcoliamo P(X <= 2) # utilizzando "pnorm" pnorm(2,mu,sigma) # ma e0 ovvio che sia 1/2, x=2 non è solo la madia ma anche la mediana!!! # calcoliamo P(X > 4)...posso pensare al complemento ad uno dell'evento negato... 1-pnorm(4,mu,sigma) # oppure pnorm(4,mu,sigma,lower.tail=F) # P(X >= -1) che in effetti è pari a P(X > -1) # infatti, essendo la v.a. continua, la presenza o meno dell'uguale NON FA DIFFERENZA! 1 - pnorm(-1,mu,sigma) # oppure pnorm(-1,mu,sigma,lower.tail = FALSE) # Disegniamo la funzione di ripartizione per diverse combinazioni di parametri # e sovrapponiamole nello stesso grafico curve(pnorm(x,0,2),xlim=c(-10,10),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="F.r di una v.a. normale \n al variare del parametro mu (sigma=2)") curve(pnorm(x,2,2),col="blue",add=T,lwd=2) curve(pnorm(x,4,2),col="red",add=T,lwd=2) curve(pnorm(x,-2,2),col="green",add=T,lwd=2) curve(pnorm(x,-4,2),col="violet",add=T,lwd=2) legend(-10,1,legend=paste("mu",c(0,2,4,-2,-4),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # Per la funz. di rip. non si hanno delle funzioni lisce proprio in virtù # dell'assoluta continuità della v.a. curve(pnorm(x,0,1),xlim=c(-12,12),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="F.r di una v.a. normale \n al variare del parametro sigma (mu=0)") curve(pnorm(x,0,2),col="blue",add=T,lwd=2) curve(pnorm(x,0,3),col="red",add=T,lwd=2) curve(pnorm(x,0,4),col="green",add=T,lwd=2) curve(pnorm(x,0,5),col="violet",add=T,lwd=2) legend(-10,1,legend=paste("sigma",1:5,sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione N(mu,sigma^2) mu=3 sigma=2 qq=seq(0.25,0.75,length=3) qnorm(qq,mu,sigma) # simulare 1000 valori da N(mu,sigma) sim.norm=rnorm(1000,mu,sigma) par(mfrow=c(2,1)) # Costruiamo l'istogramma per il campione simulato e sovrapponiamo # il modello teorico normale (o popolazione) da cui sono state effettuate # le estrazioni # Ricordiamo che l'opzione prob=T e' essenziale a questo scopo hist(sim.norm,50,prob=T,main="numero simulazioni = 1000",xlim=c(-5,10)) curve(dnorm(x,mu,sigma),col="red",add=T) # aumentiamo le simulazioni... sim.norm=rnorm(100000,mu,sigma) hist(sim.norm,100,prob=T,main="numero simulazioni = 100000",xlim=c(-5,10)) curve(dnorm(x,mu,sigma),col="red",add=T) # verifichiamo empiricamente che la media sia mu mu=3 # mu=3--> E(X)=3 mean(sim.norm) #----------------------------------------------------------------- # * variabile aleatoria chi quadrato con "k" gdl #----------------------------------------------------------------- # X v.a. continua con supporto [0,+inf] # E(X) = k # Var(X) = 2*k k=10 # funzione di densita' nel punto x=3 x=3 1/(gamma(k/2)*2^(k/2))*x^(k/2-1)*exp(-x/2) # oppure attraverso la funzione dchisq() dchisq(3,df=k) # densità al variare di k curve(dchisq(x,5),xlim=c(0,60),ylim=c(0,0.2),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. chi-quadro \n al variare dei g.d.l.") curve(dchisq(x,10),col="blue",add=T,lwd=2) curve(dchisq(x,15),col="red",add=T,lwd=2) curve(dchisq(x,20),col="green",add=T,lwd=2) curve(dchisq(x,30),col="violet",add=T,lwd=2) legend(5,0.4,legend=paste("k",c(5,10,15,20,30),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # FUNZIONE DI RIPARTIZIONE # P(X <= 15) # utilizzando "pchisq" k=20 pchisq(15,k) # P(X > 30) 1 - pchisq(30,k) # oppure pchisq(30,k,lower.tail = FALSE) # Rappresentiamo la f.r. per diversi valori dei g.d.l. curve(pchisq(x,5),xlim=c(0,60),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. chi-quadro \n al variare dei g.d.l.") curve(pchisq(x,10),col="blue",add=T,lwd=2) curve(pchisq(x,15),col="red",add=T,lwd=2) curve(pchisq(x,20),col="green",add=T,lwd=2) curve(pchisq(x,30),col="violet",add=T,lwd=2) legend(45,0.8,legend=paste("k",c(5,10,15,20,30),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione Chisq(k) k=20 qq=seq(0.25,0.75,length=3) qchisq(qq,k) # simulare 1000 valori da chisq(k) sim.chisq=rchisq(1000,k) par(mfrow=c(2,1)) hist(sim.chisq,50,prob=T,ylim=c(0,0.08),main="numero simulazioni = 1000",xlim=c(0,40)) curve(dchisq(x,k),col="red",add=T) # aumentiamo le simulazioni... sim.chisq=rchisq(100000,k) hist(sim.chisq,50,prob=T,ylim=c(0,0.08),main="numero simulazioni = 100000",xlim=c(0,40)) curve(dchisq(x,k),col="red",add=T) # verifichiamo empiricamente che la media sia k # k=20--> E(X)=20 k=20 mean(rchisq(100000,k)) # k --> E(X)=50 k=50 mean(rchisq(100000,k)) # Approssimazione normale: chisq(k)-->N(k,2k) # quando k è elevato #k=30 curve(dchisq(x,30),xlim=c(0,60),ylim=c(0,0.2),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. chi-quadro \n al variare dei g.d.l.") curve(dnorm(x,30,sqrt(2*30)),col="blue",add=T,lwd=2) # k=60 curve(dchisq(x,60),xlim=c(0,120),ylim=c(0,0.04),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. chi-quadro \n al variare dei g.d.l.") curve(dnorm(x,60,sqrt(2*60)),col="blue",add=T,lwd=2) # k=200 curve(dchisq(x,200),xlim=c(0,400),ylim=c(0,0.02),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. chi-quadro \n al variare dei g.d.l.") curve(dnorm(x,200,sqrt(2*200)),col="blue",add=T,lwd=2) #----------------------------------------------------------------- # * variabile aleatoria T di Student con "k" gdl #----------------------------------------------------------------- # X v.a. continua con supporto R # E(X) = 0 # Var(X) = k/(k-2) definita per k>2 k=10 # funzione di densita' nel punto 0.2 x=0.2 gamma((k+1)/2)/gamma(k/2)/((k*pi)^(1/2))/(((1+(x^2)/k))^((k+1)/2)) # oppure attraverso la funzione dt() dt(3,df=k) # funzione di densita' nel punto -30 dt(-30,df=k) # molto piccolo ma in effetti stiamo sulla coda sinistra! # densita' al variare di k curve(dt(x,2),xlim=c(-5,5),ylim=c(0,0.5),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. t di student \n al variare dei g.d.l.") curve(dt(x,5),col="blue",add=T,lwd=2) curve(dt(x,20),col="red",add=T,lwd=2) legend(2.5,0.4,legend=c("k=2","k=5","k=20"),lty=1,lwd=2,col=c("black","blue","red")) # FUNZIONE DI RIPARTIZIONE # P(X <= 15) # utilizzando "pt" k=20 pt(15,k) # P(X > 3) 1 - pt(3,k) # oppure pt(3,k,lower.tail = FALSE) curve(pt(x,2),xlim=c(-5,5),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. t di student \n al variare dei g.d.l.") curve(pt(x,5),col="blue",add=T,lwd=2) curve(pt(x,20),col="red",add=T,lwd=2) legend(-4,1,legend=c("k=2","k=5","k=20"),lty=1,lwd=2,col=c("black","blue","red")) # quartili della distribuzione t(k) k=20 qq=seq(0.25,0.75,length=3) qt(qq,k) # simulare 1000 valori da t(k) sim.t=rt(1000,k) # confrontiamo con la densità teorica par(mfrow=c(2,1)) hist(sim.t,50,prob=T,ylim=c(0,0.6),main="numero simulazioni = 1000",xlim=c(-4,4)) curve(dt(x,k),col="red",add=T) # aumentiamo le simulazioni... sim.t=rt(100000,k) hist(sim.t,50,prob=T,ylim=c(0,0.6),main="numero simulazioni = 100000",xlim=c(-4,4)) curve(dt(x,k),col="red",add=T) # Vediamo un modo alternativo di simulare da una t di student con k g.d.l. # tenendo conto che essa e' pari al rapporto tra una normale standardizzata e # la radice di un chi quadro diviso per i suoi k g.d.l., ove la normale e la # chi quadro sono tra loro indipendenti # Dunque, prima simulo dalla normale std camp.norm=rnorm(10000) # poi simulo da una chi quadro con opportuni g.d.l. camp.chi=rchisq(10000,k) # il campione dalla t lo ottengo col seguente rapporto camp.t=camp.norm/(sqrt(camp.chi/k)) hist(camp.t,30,prob=T,ylim=c(0,0.6),main="numero simulazioni = 10000",xlim=c(-4,4)) curve(dt(x,k),col="red",add=T) # Approssimazione Normale(0,1) per k --> +inf # Lo possiamo verificare confrontando le densita' ... curve(dnorm(x),xlim=c(-5,5),ylim=c(0,0.5),lwd=2,lty=2,xlab="x",ylab="f(x)",main="Approssimazione normale della t di Student") curve(dt(x,2),col="blue",add=T,lwd=2) curve(dt(x,20),col="red",add=T,lwd=2) legend(3,0.5,legend=c("N(0,1)","t(2)","t(20)"),lty=c(2,1,1),lwd=2,col=c("black","blue","red")) # ...o in modo equivalente le f.r. curve(pnorm(x),xlim=c(-5,5),ylim=c(0,1),lwd=2,lty=2,xlab="x",ylab="P(X<=x)",main="Approssimazione normale della t di Student") curve(pt(x,2),col="blue",add=T,lwd=2) curve(pt(x,20),col="red",add=T,lwd=2) legend(-4,1,legend=c("N(0,1)","t(2)","t(20)"),lty=c(2,1,1),lwd=2,col=c("black","blue","red")) #----------------------------------------------------------------- # * variabile aleatoria uniforme nell'intervallo [a,b] #----------------------------------------------------------------- # X v.a. continua tale per cui la f.d. e' costante in tutto l'intervallo # E(X)=(a+b)/2 # Var(X)=(b-a)^2/12 # di default v.a. uniforme in [0,1] a=-1 b=1 f.d.in.0=1/(b-a) * (a<=0 & 0<=b) f.d.in.0 f.d.in.0.5=1/(b-a) * (a<=0.5 & 0.5<=b) f.d.in.0.5 f.d.in.2=1/(b-a) * (a<=2 & 2<=b) f.d.in.2 # piu' facilmente utilizzando la funzione dunif dunif(x = 0, min = a, max = b) dunif(.5,a,b) dunif(2,a,b) # altro esempio a=0 b=1/2 f.d.in.0=1/(b-a) * (a<=0 & 0<=b) # funzione di densita' >0 (possibile anche che sia maggiore di 1...) f.d.in.0 f.d.in.0.5=1/(b-a) * (a<=0.5 & 0.5<=b) f.d.in.0.5 f.d.in.2=1/(b-a) * (a<=2 & 2<=b) f.d.in.2 # piu' facilmente utilizzando la funzione dunif dunif(0,a,b) dunif(0.5,a,b) dunif(2,a,b) # densita' al variare dei parametri curve(dunif(x,0,1),xlim=c(-5,5),ylim=c(0,2),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. uniforme \n al variare dei parametri a e b") curve(dunif(x,0,2),col="blue",add=T,lwd=2) curve(dunif(x,-1,0),col="red",add=T,lwd=2) legend(1,2,legend=c("a=0 , b=1","a=0 , b=2","a=-1 , b=0"),lty=1,lwd=2,col=c("black","blue","red")) # FUNZIONE DI RIPARTIZIONE # P(X <= 0.3) per uniforme in [0,1] # utilizzando "punif" punif(0.3) # a=0 ; b=1 # X~Unif(0.1) --> P(X 0.4) per uniforme in [-1,1] 1 - punif(0.4,a,b) # oppure punif(0.4,a,b,lower.tail = FALSE) curve(punif(x,0,1),xlim=c(-2,3),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. uniforme \n al variare dei parametri a e b") curve(punif(x,0,2),col="blue",add=T,lwd=2) curve(punif(x,-1,0),col="red",add=T,lwd=2) legend(1,2,legend=c("a=0 , b=1","a=0 , b=2","a=-1 , b=0"),lty=1,lwd=2,col=c("black","blue","red")) # quartili della distribuzione uniforme [a,b] qq=seq(0.25,0.75,length=3) qunif(qq) # di default a=0 e b=1 qunif(qq,a,b) # simulare 1000 valori da unif(a,b) a=-1 b=1 sim.unif=runif(1000,a,b) par(mfrow=c(2,1)) hist(sim.unif,50,prob=T,main="numero simulazioni = 1000") curve(dunif(x,a,b),col="red",add=T) # aumentiamo le simulazioni... sim.unif=runif(100000,a,b) hist(sim.unif,50,prob=T,main="numero simulazioni = 100000") curve(dunif(x,a,b),col="red",add=T) # verifichiamo empiricamente che la media sia (a+b)/2 # a=0 ; b=1 --> E(X)=0.5 mean(runif(100000)) # a=-1 ; b=1 --> E(X)=0 mean(runif(100000,a,b)) #----------------------------------------------------------------- # * variabile aleatoria beta di parametri "a" e "b" #----------------------------------------------------------------- # X v.a. continua con supporto [0,1] # a,b > 0 # E(X) = a/(a+b) # Var(X) = ab/((a+b+1)*(a+b)^2) a=2 b=2 # funzione di densita' nel punto 0.2 (0.2^(a-1)*(1-0.2)^(b-1)) / beta(a,b) # dove beta(a,b) e' la funzione beta da non confondersi con la distribuzione beta! # oppure attraverso la funzione dbeta() dbeta(0.2,shape1=a,shape2=b) # NOTA: a,b=1 --> uniforme(0,1) dbeta(0.2,1,1) dbeta(0.3,1,1) dunif(0.2) # f.d. costante e pari a 1 dunif(0.3) # vediamo come varia la densita' beta al variare dei parametri a e b curve(dbeta(x,1,1),xlim=c(0,1),ylim=c(0,4),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. beta \n al variare dei parametri a e b") curve(dbeta(x,2,2),col="blue",add=T,lwd=2) curve(dbeta(x,0.5,0.5),col="red",add=T,lwd=2) curve(dbeta(x,0.5,2),col="green",add=T,lwd=2) curve(dbeta(x,2,0.5),col="violet",add=T,lwd=2) legend(0.6,4,legend=c("a=1 , b=1","a=2 , b=2","a=0.5 , b=0.5","a=0.5 , b=2","a=2 , b=0.5"),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # a,b > 1 curve(dbeta(x,2,2),xlim=c(0,1),ylim=c(0,4),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. beta \n al variare dei parametri a e b (a,b>1)") curve(dbeta(x,5,5),col="blue",add=T,lwd=2) curve(dbeta(x,2,5),col="red",add=T,lwd=2) curve(dbeta(x,5,2),col="green",add=T,lwd=2) legend(0.6,4,legend=c("a=2 , b=2","a=5 , b=5","a=2 , b=5","a=5 , b=2"),lty=1,lwd=2,col=c("black","blue","red","green")) #a,b<1 curve(dbeta(x,0.5,0.5),xlim=c(0,1),ylim=c(0,4),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. beta \n al variare dei parametri a e b (a,b>1)") curve(dbeta(x,0.1,0.1),col="blue",add=T,lwd=2) curve(dbeta(x,0.5,0.1),col="red",add=T,lwd=2) curve(dbeta(x,0.1,0.5),col="green",add=T,lwd=2) legend(0.6,4,legend=c("a=0.5 , b=0.5","a=0.1 , b=0.1","a=0.5 , b=0.1","a=0.1 , b=0.5"),lty=1,lwd=2,col=c("black","blue","red","green")) # FUNZIONE DI RIPARTIZIONE # P(X <= 0.3) # utilizzando "pbeta" pbeta(0.3,a,b) # P(X > 0.4) 1 - pbeta(0.4,a,b) # oppure pbeta(0.4,a,b,lower.tail = FALSE) curve(pbeta(x,1,1),xlim=c(0,1),ylim=c(0,1),lwd=2,xlab="x",ylab="P(X<=x)",main="Funzione di ripartizione di una v.a. beta \n al variare dei parametri a e b") curve(pbeta(x,2,2),col="blue",add=T,lwd=2) curve(pbeta(x,0.5,0.5),col="red",add=T,lwd=2) curve(pbeta(x,0.5,2),col="green",add=T,lwd=2) curve(pbeta(x,2,0.5),col="violet",add=T,lwd=2) legend(-0.01,1,legend=c("a=1 , b=1","a=2 , b=2","a=0.5 , b=0.5","a=0.5 , b=2","a=2 , b=0.5"),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione beta(a,b) a=2 b=2 qq=seq(0.25,0.75,length=3) qbeta(qq,a,b) # simulare 1000 valori da beta(a,b) sim.beta=rbeta(1000,a,b) par(mfrow=c(2,1)) hist(sim.beta,50,prob=T,main="numero simulazioni = 1000") curve(dbeta(x,a,b),col="red",add=T) # aumentiamo le simulazioni... sim.beta=rbeta(100000,a,b) hist(sim.beta,50,prob=T,main="numero simulazioni = 100000") curve(dbeta(x,a,b),col="red",add=T) # verifichiamo empiricamente che la media sia a/(a+b) # a=2 ; b=2 --> E(X)=0.5 mean(rbeta(100000,2,2)) # a=3 ; b=6 --> E(X)=1/3 a=3 b=6 mean(rbeta(100000,a,b)) # a=6 ; b=3 --> E(X)=2/3 a=6 b=3 mean(rbeta(100000,a,b)) #----------------------------------------------------------------- # * variabile aleatoria F di Fisher con "n" e "m" gdl #----------------------------------------------------------------- # X v.a. continua con supporto [0,+inf] # n,m > 0 # E(X) = m/(m-2) # funzione di densita' nel punto x=2 n=15 m=10 x=2 sqrt((((n*x)^n)*m^m)/((n*x+m)^(n+m)))/(x*beta((n/2),(m/2))) # oppure attraverso la funzione dt() df(x,df1=n,df2=m) curve(df(x,2,2),xlim=c(0,10),ylim=c(0,1),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. F di Fisher \n al variare dei g.d.l.") curve(df(x,5,5),col="blue",add=T,lwd=2) curve(df(x,20,20),col="red",add=T,lwd=2) curve(df(x,3,50),col="orange",add=T,lwd=2) curve(df(x,50,3),col="green",add=T,lwd=2) legend(7,1,legend=c("n=m=2","n=m=5","n=m=20","n=3 ; m=50","n=50 ; m=3"),lty=1,lwd=2,col=c("black","blue","red","orange","green")) # FUNZIONE DI RIPARTIZIONE # P(X <= 2) # utilizzando "pf" n=20 m=15 pf(2,n,m) # P(X > 3) 1 - pf(3,n,m) # oppure pf(3,n,m,lower.tail = FALSE) curve(pf(x,2,2),xlim=c(0,10),ylim=c(0,1),lwd=2,xlab="x",ylab="f(x)",main="Funzione di ripartizione di una v.a. F di Fisher \n al variare dei g.d.l.") curve(pf(x,5,5),col="blue",add=T,lwd=2) curve(pf(x,20,20),col="red",add=T,lwd=2) curve(pf(x,3,50),col="orange",add=T,lwd=2) curve(pf(x,50,3),col="green",add=T,lwd=2) legend(7,1,legend=c("n=m=2","n=m=5","n=m=20","n=3 ; m=50","n=50 ; m=3") # quartili della distribuzione f(k) n=20 m=15 qq=seq(0.25,0.75,length=3) qf(qq,n,m) # simulare 1000 valori da f(k) sim.f=rf(1000,n,m) par(mfrow=c(2,1)) hist(sim.f,100,prob=T,ylim=c(0,1.5),main="numero simulazioni = 1000",xlim=c(0,6)) curve(df(x,n,m),col="red",add=T) # aumentiamo le simulazioni... sim.f=rf(100000,n,m) hist(sim.f,100,prob=T,ylim=c(0,1.5),main="numero simulazioni = 100000",xlim=c(0,6)) curve(df(x,n,m),col="red",add=T) # verifichiamo empiricamente che la media sia m/(m-2) # m=20--> E(X)=m/(m-2)=1.111 n=100 m=20 mean(rf(100000,n,m)) # m=5 --> E(X)=1.667 n=100 m=5 mean(rf(100000,n,m)) #----------------------------------------------------------------- # * variabile aleatoria esponenziale di parametro "mu" #----------------------------------------------------------------- # X v.a. continua con supporto [0,+inf] # mu > 0 # E(X) = 1/mu # Var(X) = 1/mu^2 mu=2 # funzione di densita' nel punto 0.2 mu*exp(-mu*0.2) # oppure attraverso la funzione dexp() dexp(0.2,rate=mu) # vediamo come varia la densita' espenenziale al variare del parametro mu curve(dexp(x,1),xlim=c(0,8),ylim=c(0,3),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. esponenziale \n al variare del parametro mu") curve(dexp(x,2),col="blue",add=T,lwd=2) curve(dexp(x,4),col="red",add=T,lwd=2) curve(dexp(x,0.5),col="green",add=T,lwd=2) curve(dexp(x,0.7),col="violet",add=T,lwd=2) legend(5,3,legend=paste("mu",c(1,2,4,0.5,0.7),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # FUNZIONE DI RIPARTIZIONE # P(X <= 1) # utilizzando "pexp" mu=2 pexp(1,mu) # P(X > 0.5) 1 - pexp(0.5,mu) # oppure pexp(0.5,mu,lower.tail = FALSE) curve(pexp(x,1),xlim=c(0,8),ylim=c(0,1),lwd=2,xlab="x",ylab="f(x)",main="Funzione di ripartizione di una v.a. esponenziale \n al variare del parametro mu") curve(pexp(x,2),col="blue",add=T,lwd=2) curve(pexp(x,4),col="red",add=T,lwd=2) curve(pexp(x,0.5),col="green",add=T,lwd=2) curve(pexp(x,0.7),col="violet",add=T,lwd=2) legend(6,0.85,legend=paste("mu",c(1,2,4,0.5,0.7),sep="="),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione exp(mu) qq=seq(0.25,0.75,length=3) qexp(qq,mu) # simulare 1000 valori da exp(mu) mu=5 sim.exp=rexp(1000,mu) par(mfrow=c(2,1)) hist(sim.exp,prob=T,ylim=c(0,6),50,main="numero simulazioni = 1000",xlim=c(0,6)) curve(dexp(x,mu),col="red",add=T) # aumentiamo le simulazioni... sim.exp=rexp(100000,mu) hist(sim.exp,prob=T,ylim=c(0,6),50,main="numero simulazioni = 100000",xlim=c(0,6)) curve(dexp(x,mu),col="red",add=T) # verifichiamo empiricamente che la media sia 1/mu # mu=2 --> E(X)=0.5 mu=2 mean(rexp(100000,mu)) # mu=0.5 --> E(X)=2 mu=0.5 mean(rexp(100000,mu)) #----------------------------------------------------------------- # * variabile aleatoria gamma di parametri "alpha" e "beta" #----------------------------------------------------------------- # X v.a. continua con supporto [0,+inf] # alpha,beta > 0 # E(X) = alpha/beta # Var(X) = alpha/beta^2 alpha=2 beta=3 # funzione di densita' nel punto 0.5 x=0.5 beta^alpha/gamma(alpha)*x^(alpha-1)*exp(-beta*x) # dove gamma(a,b) e' la funzione gamma da non confondersi # con la densita' gamma! # oppure attraverso la funzione dgamma() dgamma(0.5,shape=alpha,rate=beta) # densita' al variare dei parametri curve(dgamma(x,1,1),xlim=c(0,10),ylim=c(0,2),lwd=2,xlab="x",ylab="f(x)",main="Densita' di una v.a. esponenziale \n al variare dei parametri") curve(dgamma(x,2,2),col="blue",add=T,lwd=2) curve(dgamma(x,10,10),col="red",add=T,lwd=2) curve(dgamma(x,5,1),col="green",add=T,lwd=2) curve(dgamma(x,1,5),col="violet",add=T,lwd=2) legend(5,1.5,legend=c("alpha=1 , beta=1","alpha=2 , beta=2","alpha=10 , beta=10","alpha=5 , beta=1","alpha=1 , beta=5"),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # FUNZIONE DI RIPARTIZIONE # P(X <= 1) # utilizzando "pgamma" alpha=2 beta=3 pgamma(1,alpha,beta) # P(X > 0.5) 1 - pgamma(0.5,alpha,beta) # oppure pgamma(0.5,alpha,beta,lower.tail = FALSE) curve(pgamma(x,1,1),xlim=c(0,10),ylim=c(0,1),lwd=2,xlab="x",ylab="f(x)",main="Funzione di ripartizione di una v.a. esponenziale \n al variare dei parametri") curve(pgamma(x,2,2),col="blue",add=T,lwd=2) curve(pgamma(x,10,10),col="red",add=T,lwd=2) curve(pgamma(x,5,1),col="green",add=T,lwd=2) curve(pgamma(x,1,5),col="violet",add=T,lwd=2) legend(6,0.4,legend=c("alpha=1 , beta=1","alpha=2 , beta=2","alpha=10 , beta=10","alpha=5 , beta=1","alpha=1 , beta=5"),lty=1,lwd=2,col=c("black","blue","red","green","violet")) # quartili della distribuzione exp(mu) qq=seq(0.25,0.75,length=3) qgamma(qq,alpha,beta) # simulare 1000 valori da exp(mu) alpha=2 beta=3 sim.gamma=rgamma(1000,alpha,beta) par(mfrow=c(2,1)) hist(sim.gamma,100,prob=T,ylim=c(0,1.5),xlim=c(0,3),main="numero simulazioni = 1000") curve(dgamma(x,alpha,beta),col="red",add=T) # aumentiamo le simulazioni... sim.gamma=rgamma(100000,alpha,beta) hist(sim.gamma,100,prob=T,ylim=c(0,1.5),xlim=c(0,3),main="numero simulazioni = 100000") curve(dgamma(x,alpha,beta),col="red",add=T) # verifichiamo empiricamente che la media sia alpha/beta # alpha=2 beta=3 --> E(X)=2/3 alpha=2 beta=3 mean(rgamma(100000,alpha,beta)) # alpha=5 beta=3 --> E(X)=5/3 alpha=5 beta=3 mean(rgamma(100000,alpha,beta)) # v.a. chi-quadro come caso particolare della v.a. gamma # gamma(alpha,beta) t.c. alpha=k/2 & beta=1/2 ---> chiq(k) k=10 hist(rchisq(100000,k),100,prob=T,main="Chi-quadro come caso particolare della Gamma") curve(dgamma(x,k/2,1/2),lwd=2,col="red",xlim=c(0,30),add=T) # v.a. esponenziale come caso particolare della v.a. gamma # gamma(alpha,beta) t.c. alpha=1 ---> exp(beta) beta=10 hist(rexp(100000,beta),100,prob=T,main="Esponenziale come caso particolare della Gamma") curve(dgamma(x,1,beta),lwd=2,col="red",xlim=c(0,3),add=T)