############################### # Esercitazione 4: 31/10/2024 # # Laboratorio in R # ############################### getwd() # Vedere l'attuale directory di lavoro # Modificare la directory di lavoro se volete caricare o salvare # qualcosa in questa specifica cartella di lavoro setwd("!!!!! !!!!") ls() # Vedere gli oggetti presenti nell'ambiente rm(list=ls()) # Per pulire l'ambiente ###################################################################### ## ALCUNI COMANDI UTILI per lavorare con scalari, vettori e matrici ## ###################################################################### ########### # Scalari # ########### a <- 2 b <- 6 # Somma a + b somma <- a + b somma ls() # Sottrazione a - b # Moltiplicazione a * b # Divisione a/b # Elevamento a potenza a^3 # Operatori di confronto a > 3 a < 3 # Operatori logici # intersezione (a > 1) & (b < 7) (a > 2) & (b < 7) (a > 2) & (b < 6) # unione (a > 2) | (b < 7) # negazione !(a > 2) # Funzioni matematiche log(a) # logaritmo naturale exp(a) # esponenziale sqrt(a) # radice quadrata sin(a) # seno cos(a) # coseno abs(-a) # valore assoluto min(a,b) # minimo max(a,b) # massimo ########### # Vettori # ########### v <- c(100, 18, 41, 21, 53) v # Accesso agli elementi del vettore v[2] # Secondo elemento del vettore v[1:3] # Primi tre elementi del vettore v > 50 v[v > 20 & v < 80] v * 3 length(v) # Funzioni matematiche sqrt(v) log(v) min(v) max(v) sum(v) cumsum(v) help("mean") #per vedere la pagina di aiuto con la descrizione della funzione mean mean(v) sortvett <- sort(v) sortvett median(v) # Operazioni tra vettori v2 <- c(1, 3, 5, 7, 9) v + v2 v * v2 v / v2 v%*%v2 ########### # Matrici # ########### A <- matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3, byrow = TRUE) A A[1, 2] A[1,] A[, 2] A[, 2 : 3] t(A) sum(A) max(A) min(A) 3 * A A/3 B <- matrix(c(2,4,6,8,10,12), nrow = 2, ncol = 3, byrow = T) A + B A - B A * B A / B A%*%B # da errore ovviamente A%*%t(B) ########################################################################## # Semplice esempio di funzione per rappresentare i dati: funzione plot() # ########################################################################## x <- c(0.5, 0.5, 2, 4, 6, 0.2, 5, 4.1) y <- c(0.5, 1.5, 3.2, 8, 6, 1.2, 3, 3.5) plot(x, y, pch = 21, bg = "red") ###################################################### # Variabili aleatorie e distribuzioni di probabilità # ###################################################### ################################ # Variabili aleatorie discrete # ################################ ## Distribuzione binomiale X ~ Bin(n = 20, p = 0.5) # Valori del supporto: sequenza di interi da 0 a 20 di passo 1 x <- seq(from = 0, to = 20, by = 1) x # Ugualmente x <- 0 : 20 x # Otteniamo P(X=x) a mano # (nota: choose permette di ottenere il coefficiente binomiale) pbin <- choose(20, x) * 0.5^x * (1 - 0.5)^(20 - x) pbin # Verifichiamo che sommi a 1 sum(pbin) # Calcoliamo la probabilità P(X=5) choose(20, 5) * 0.5^5 * (1 - 0.5)^(20 - 5) # Otteniamo il grafico della funzione di probabilità plot(x, pbin, type = "h", main = "Binomiale(20, 0.5)", xlab = "x", ylab = "p(x)") # Costruiamo una funzione che ci permette di ottenere le probabilità # passando in argomento il numero di prove e la probabilità di successo f.bin <- function(x, n, p){ return(choose(n, x) * p^x * (1 - p)^(n - x)) } # Confronto tra due distribuzioni binomiali x <- 0 : 20 y0 <- f.bin(x, n = 20, p = 0.3) y1 <- f.bin(x, n = 20, p = 0.6) plot(x, y0, type = "h", ylim = c(0, 0.2), xlab = "x", ylab = "p(x)") points(x - 0.1, y1, type = "h", col = 2) # Utilizzo di funzioni predefinite di R # P(X=2) con X ~ Bin(n = 6, p = 0.55) dbinom(2, 6, 0.55) f.bin(2, 6, 0.55) # P(X>2) con X ~ Bin(n = 6, p = 0.55) 1 - pbinom(2, 6, 0.55) 1 - sum(dbinom(0 : 2, 6, 0.55)) pbinom(2, 6, 0.55, lower.tail = FALSE) # Confrontare le funzioni di 2 binomiali per n fiissato # e diversi valori di p par(mfrow = c(1, 2)) x <- 0 : 20 n <- 20 plot(x, dbinom(x, n, 0.1), type = "h", xlab = "x", ylab = "P(x)", ylim = c(0, 0.35), main = "Bin(n = 20, p = 0.1)") plot(x, dbinom(x, n, 0.3),type = "h", xlab = "x", ylab = "P(x)", ylim = c(0, 0.35), main = "Bin(n = 20, p = 0.3)") # Approssimazione Binomiale - Poisson (sull'esercizio 5 dell'esercitazione 3) n <- 2000 p <- 0.03 x <- seq(30, 90, by = 1) # Si produce il grafico solo per le x comprese tra 30 e 90 par(mfrow = c(1, 1)) plot(x, dpois(x, n * p), type = "h", ylab =" ") points(x + 0.2, dbinom(x, n, p), type = "h", col = "red") # P(T <=50) T ~ Bin(n = 2000, p = 0.03) pbinom(50, n, p) # usando l'approssimazione con la Poisson ppois(50, n* p) #### Generazione di numeri casuali ## Da una binomiale y <- rbinom(100, size = 1, prob = 0.5) sum(y)/100 mean(y) ## Da una Poisson cpois1 <- rpois(10, 1) # media=1 mean(cpois1) cpois2 <- rpois(10, 2) # media=2 mean(cpois2) cpois3 <- rpois(10, 20) # media=20 mean(cpois3) ############################## # Risultati diversi: perchè? # vpois <- rpois(100, 3) mean(vpois) vpois <- rpois(100, 3) mean(vpois) # Fissare un seme set.seed(13) vpois <- rpois(100, 3) mean(vpois) set.seed(13) vpois <- rpois(100, 3) mean(vpois) # Generiamo 500 valori casuali da una distribuzione di Poisson con parametro 5 n <- 500 set.seed(13) dati.pois <- rpois(n, 5) #500 valori dalla Po(5) # Visualizziamo la tabelle delle frequenze tab <- table(dati.pois) tab # Verifichiamo l'adattamento mediante bar plot plot(tab/n) m <- max(dati.pois) x <- 0 : m points(x + 0.1, dpois(x, 5), type = "h", col = 2) ################################ # Variabili aleatorie continue # ################################ ## Distribuzione Normale # Funzione di densità valutata in 1 per N(0,1) e N(0, 1.5) dnorm(1) dnorm(1, mean = 1, sd = 1.5) # Funzione di ripartizione valutata in 1 per N(0,1) e N(0, 1.5) pnorm(1) pnorm(1, mean = 1, sd = 1.5) # Grafico funzioni di densità x <- seq(-6, 10, length = 100) plot(x, dnorm(x), type = "l", xlab = "x", ylab = "densità") curve(dnorm(x, mean = 1, sd = 1.5), add = TRUE, lty = 2, col = 2) points(x, dnorm(x, mean = 3, sd = 2), type = "l", col = 3, lwd = 2) # quantile di ordine 0.10 qnorm(0.1) -qnorm(0.9) # Grafico funzioni di ripartizione sigma <- 1.3 x <- seq(-10, 10, length = 100) plot(x, pnorm(x, mean = -1, sd = sigma), type = "l", xlab = "x", ylab = "F(x)") curve(pnorm(x, mean = 0, sd = sigma), add = TRUE, col = 2) points(x, pnorm(x, mean = 3, sd = sigma), type = "l", col = 3) points(x, pnorm(x, mean = 3, sd = sigma * 2), type = "l", col = 4) # Generazione di campioni casuali da distribuzione N(10,25) e verifica adattamento set.seed(5) x <- rnorm(800, 10, 5) #simuliamo da una N(10, 25) hist(x, freq = FALSE, breaks = 30, xlab = " ", ylab = "densità", main = "") curve(dnorm(x, 10, 5), add = T, col = 4, lwd = 2) # Funzione di ripartizione empirica set.seed(9) Z <- rnorm(150) #N(0,1) edf_norm <- ecdf(Z) # f.d.r empirica edf_norm # Grafico f.d.r. empirica e teorica plot(edf_norm, verticals = TRUE, do.p=FALSE) curve(pnorm(x), add = T, col = 2) # Sull'esempio delle assicurazioni automobilistiche dove ci chiedevamo n <- 2000 p <- 0.03 pnorm(50, mean = n * p, sd = sqrt(n * p * (1 - p)) ) #################################### # Simulazione e metodi Monte Carlo # #################################### # Illustrazione di un ciclo for #numero iterazioni nsim <- 5 # in 'out' è memorizzato l'output di ogni simulazione out <- vector(mode = "numeric", length = nsim) for(i in 1 : nsim){ # Qui dichiaro le istruzioni da eseguire ad ogni iterazione #stampa il valore corrente dell'iterazione cat("iter", i, "") #il valore corrente viene memorizzato nella posizione i-ma out[i] = i } out ################################ # Somma di variabili aleatorie # ################################ ############################################# ## Somme di v.a. esponenziali indipendenti ## ############################################# set.seed(2) nsim <- 1000 #numero di iterazioni n <- 10 lambda <- 0.1 #x conterrà la somma dei valori estratti x <- vector(mode = "numeric", length = nsim) for(i in 1 : nsim){ campione <- rexp(n, lambda) x[i] <- sum(campione) } hist(x, prob = TRUE, ylim = c(0, 0.015), ylab = "densità", main = "Confronto valori simulati e teorici", breaks = 30, cex.main = 0.9) curve(dgamma(x, n, lambda), type = "l", col = 2, add = T, lwd = 2) ########################################## ## Somma di variabili aleatorie normali ## ########################################## set.seed(1) nsim <- 10000 n <- 10 mu <- 2 sigma <- 1 S <- vector(mode = "numeric", length = nsim) for (i in 1 : nsim){ S[i] <- sum(rnorm(n, mean = mu, sd = sigma)) } hist(S, prob = TRUE, breaks = 30, ylim = c(0, 0.18)) curve(dnorm(x, mean = n * mu, sd = sqrt(n)), type = "l", add = T, col = 3, lwd = 2) #################################################### ## TLC: Somma di variabili aleatorie rettangolari ## #################################################### set.seed(100) nsim <- 1000 n <- 20 R <- vector(mode = "numeric", length = nsim) for(i in 1 : nsim){ R[i] <- sum(runif(n)) } str(R) hist(R, prob = T, breaks = 40, xlab = " ", ylab = " ", main = "Distribuzione della somma di R(0,1)", cex.main = .9) curve(dnorm(x, mean = n * (1/2), sd = sqrt(n/12)), type = "l", add = T, col = 3, lwd = 2) ################################################################# ## Distribuzione campionaria della media di un campione normale # ################################################################# set.seed(1) R <- 1000 # replicazioni n <- 10 mc <- vector(mode = "numeric", length = R) for(i in 1 : R){ x <- rnorm(n) mc[i] <- mean(x) } # Verifica dell'adattamento par(mfrow = c(1,3)) hist(mc, freq = FALSE, main = "", ylim = c(0,1.3)) curve(dnorm(x, mean = 0, sd = sqrt(1/n)), add = TRUE, col = 2, lwd = 2) qqplot(mc, rnorm(R, mean = 0, sd = sqrt(1/n)), ylab = " ") abline(0, 1) boxplot(mc) ################################################################################### ## Distribuzione campionaria della media di un campione non normale - rettangolare# ################################################################################### set.seed(2) R <- 1000 #numero di replicazioni n <- 30 #dimensione del campione medie <- vector(mode = "numeric", length = R) for(i in 1 : R){ x <- runif(n) medie[i] <- mean(x) } mean(medie) var(medie) par(mfrow = c(1, 1)) hist(medie, breaks = 30, prob = TRUE) curve(dnorm(x, mean = (1/2), sd = 1/sqrt(12 * n)), col = 2, lwd = 2, add = T) ############################################## ## Distribuzione campionaria della varianza ## ############################################## set.seed(33) R <- 1000 n <- 8 mu <- 5 sigma <- 2 # Che vuol dire sigma^2 = 4 S2 <- S2nc <- vector(mode="numeric", length=R) for(i in 1 : R){ x <- rnorm(n, mu, sigma) S2[i] <- var(x) S2nc[i] <- S2[i] * (n - 1)/n } # Verifichiamo l'adattamento par(mfrow=c(1,2)) hist(S2, prob = TRUE, breaks = 30, main = "S2") hist(S2nc, prob = TRUE, breaks = 30, main = "S2nc") # Confrontiamo la media dei valori generati delle due varianze campionarie mean(S2) #media del campione di S2 mean(S2nc) #media del campione di S2nc # Sfruttando il risultato teorico verifichiamo l'adattamento par(mfrow=c(1,1)) hist(S2, prob = TRUE, breaks = 30, ylab = "densità", main = " ") curve(((n - 1)/sigma^2) * dchisq(x * ((n - 1)/sigma^2), df = n - 1), add = TRUE, col = 4, lwd = 2, main = "") ############################################# ## Trasformazione integrale di probabilità ## ############################################# # Generiamo valori casuali da una distribuzione esponenziale # di parametro lambda R <- 1000 lambda <- 2 x <- rep(0, R) for(i in 1 : R) { u <- runif(1, 0, 1) x[i] <- -log(u)/lambda } hist(x, prob = TRUE, breaks = 30) curve(dexp(x, lambda), from = 1e-16, to = 4, add = T, col = "red")