###############################
# 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("!!!!! <inserire la directory> !!!!")   

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")