rm(list = ls()) #svuotare l'ambiente di lavoro
try(dev.off()) #azzerare i paramentri grafici
gc() #liberare la ram

library(ggplot2)

# 0. Comandi base ####
#ottenere la probabilità esatta in relazione al valore, P(x = X)
dexp(1,1) #esponenziale
dnorm() #normale
dbinom() #binomiale
dnbinom() #binomiale negativa
dpois() #poisson
dgamma() #gamma
dchisq() #chi-quadrato
dt() #t-student

#ottenere la probabilità cumulata in relazione al valore, P(x <= X)
# NOTA: ...lower.tail = T  se P(x > X)
pexp(1,1,lower.tail = F) #esponenziale
pnorm() #normale
pbinom() #binomiale
pnbinom() #binomiale negativa
ppois() #poisson
pgamma() #gamma
pchisq() #chi-quadrato
pt() #t-student

#ottenere il valore in relazione alla probabilità
qexp(0.6321206,1) #esponenziale
qnorm() #normale
qbinom() #binomiale
qnbinom() #binomiale negativa
qpois() #poisson
qgamma() #gamma
qchisq() #chi-quadrato
qt() #t-student

#generare campioni casuali
rexp() #esponenziale
rnorm() #normale
rbinom() #binomiale
rnbinom() #binomiale negativa
rpois() #poisson
rgamma() #gamma
rchisq() #chi-quadrato
rt() #t-student

# 1. Generazione e visualizzazione di distribuzioni ####

## 1.1. Distribuzione esponenziale ####

### 1.1.1. Generazione campione e istogramma ####
set.seed(1) #fissiamo il seme della casualità per la riproducibiltà dei risultati
x <- rexp(100) #generiamo un campione di 100 elementi
dir.create("data")
dir.create("data/campioni casuali")
save(x, file = "data/campioni casuali/esponenziale.Rdata")
plot(x) #visualizziamoli
plot(x, type = "l")
plot(x, type = "b")
plot(x, type = "o")
plot(x, type = "s")
plot(x, type = "h")
# type style (from the help):
# what type of plot should be drawn. Possible types are
# "p" for points,
# "l" for lines,
# "b" for both,
# "c" for the lines part alone of "b",
# "o" for both ‘overplotted’,
# "h" for ‘histogram’ like (or ‘high-density’) vertical lines,
# "s" for stair steps,
# "S" for other steps, see ‘Details’ below,
# "n" for no plotting.

hist(x)
hist(x, breaks = 50)
hist(x, breaks = 20)
hist(x, freq = F, breaks = 20)
#generare la distribuzione teorica
x_g <- seq(0, 5, .1)
y <- exp(-x_g) #distribuzione di probabilità
lines(x_g, y) # aggiungiamo la linea

lines(density(x), col = "blue") #guardare la density() function nell help
#non è detto sia sempre la scelta ottimale

library(ggplot2)
ggplot() +
  geom_histogram(aes(x = x, y = after_stat(density)),
                 bins = 30,
                 col = "black") +
  geom_density(aes(x), col = "black") +
  geom_line(aes(y = y, x = x_g), col = "blue")

### ESERCIZIO (1.1.1.): ####
# Creare una funzione ad un solo paramentro (numerosità del campione)
# che generi e salvi in un file esterno l'istogramma del campione

# ESEMPIO SVOLTO
# genera diverse immagini a partire aumentando la numerosità del campione
dir.create("plot/hist_exp") #creare la cartella di destinazione
hist_exp <- function(n) {
  x <- rexp(n)
  p <- ggplot() +
    geom_histogram(aes(x = x, y = after_stat(density)),
                   bins = 50,
                   col = "black") +
    geom_density(aes(x), col = "black") +
    geom_line(aes(y = y, x = x_g), col = "blue")
  pdf(paste0("plot/hist_exp/hist_exp_", n, ".pdf"))
  print(p)
  dev.off()
}
for (n in c(10 ^ 2, 10 ^ 3, 10 ^ 4, 10 ^ 5)) {
  hist_exp(n)
}

# Cosa succede se aumentiamo la numerosità del campione?


### 1.1.2. Calcolo della media, varianza, ecc. ####

load("data/campioni casuali/esponenziale.Rdata")
lambda <- 1
print(paste("media teorica:", 1 / lambda))
print(paste("media campionaria", mean(x)))
print(paste("mediana teorica:", log(2) / lambda))
print(paste("mediana campionaria:", median(x)))
print(paste("varianza teorica:", 1 / lambda ^ 2))
print(paste("varianza campionaria:", var(x)))

stat_descr_exp <- function(n, lambda) {
  x <- rexp(n, lambda)
  print("NUOVO CAMPIONE")
  print(paste("n:", n, "lambda:", lambda))
  print(paste("media teorica:", 1 / lambda))
  print(paste("media campionaria", mean(x)))
  print(paste("mediana teorica:", log(2) / lambda))
  print(paste("mediana campionaria:", median(x)))
  print(paste("varianza teorica:", 1 / lambda ^ 2))
  print(paste("varianza campionaria:", var(x)))
}

for (lambda in c(.5, 1, 3, 5, 10)) {
  for (n in c(10 ^ c(1, 3, 5))) {
    stat_descr_exp(n, lambda)
  }
}

#come creare una tabella che racchiuda i valori generati dal ciclo for
# creare un dataframe vuoto
df <- data.frame(
  n = numeric(),
  lambda = numeric(),
  media_teo = numeric(),
  media_camp = numeric(),
  mediana_teo = numeric(),
  mediana_camp = numeric(),
  var_teo = numeric(),
  var_camp = numeric()
)

for (lambda in c(.5, 1, 3, 5, 10)) {
  for (n in c(10 ^ seq(1, 5, .2))) {
    x <- rexp(n, lambda)
    df <- rbind(
      df,
      data.frame(
        n = n,
        lambda = lambda,
        media_teo = 1 / lambda,
        media_camp = mean(x),
        mediana_teo = log(2) / lambda,
        mediana_camp = median(x),
        var_teo = 1 / lambda ^ 2,
        var_camp = var(x)
      )
    )
  }
}

View(df)

ggplot(df) +
  geom_line(aes(x = n, y = media_teo, col = as.factor(lambda)), linetype =
              2) +
  geom_line(aes(x = n, y = media_camp, col = as.factor(lambda))) +
  ggtitle(label = "media teorica vs campionaria")

#ri-scalare gli assi
ggplot(df) +
  geom_line(aes(x = n, y = media_teo, col = as.factor(lambda)), linetype =
              2) +
  geom_line(aes(x = n, y = media_camp, col = as.factor(lambda))) +
  ggtitle(label = "media teorica vs campionaria") +
  scale_x_continuous(trans = "log") + ylab("")

#ri-scalare gli assi e aggiungere le etichette verticali
ggplot(df) +
  geom_line(aes(x = n, y = media_teo, col = as.factor(lambda)), linetype =
              2) +
  geom_line(aes(x = n, y = media_camp, col = as.factor(lambda))) +
  ggtitle(label = "media teorica vs campionaria") +
  scale_x_continuous(trans = "log", breaks = unique(df$n)) +
  theme(axis.text.x = element_text(
    angle = 90,
    hjust = 1,
    vjust = .5
  )) +
  scale_color_discrete(name = "lambda") + ylab("")

#ripetere per varianza e mediana
convergence_plot <- function(target_cols) {
  p <- ggplot(df) +
    geom_line(aes(
      x = n,
      y = df[, target_cols[1]],
      col = as.factor(lambda)
    ), linetype = 2) +
    geom_line(aes(
      x = n,
      y = df[, target_cols[2]],
      col = as.factor(lambda)
    )) +
    ggtitle(label = paste(names(df)[target_cols[1]], "vs campionaria")) +
    scale_x_continuous(trans = "log", breaks = unique(df$n)) +
    theme(axis.text.x = element_text(
      angle = 90,
      hjust = 1,
      vjust = .5
    )) +
    scale_color_discrete(name = "lambda") + ylab("")
  return(p)
}
for (n in c(5, 7)) {
  print(convergence_plot(c(n, n + 1)))
}

## 1.2 Distribuzione Gaussiana (per casa) ####

## 1.3 Distribuzione di Bernoulli (per casa) ####

## 1.4 Generalizzazione ad ogni distribuzione ####
hist_dist <- function(dist,n,param){
if (dist=="exp") {
  x<-rexp(n,param[1])
}
if (dist=="norm") {
  x<-rnorm(n,param[1],param[2])
}
if (dist=="bino") {
  x<-rbinom(param[1],n,param[2])
}
if (dist=="nbino") {
  x<-rnbinom(param[1],n,param[2])
}
if (dist=="pois") {
  x <- rpois(n,param[1])
}
if (dist=="gamma") {
  x <- rgamma(n,param[1])
}
if (dist=="chisq") {
  x <- rchisq(n,param[1])
}
if (dist=="t") {
  x <- rt(n,param[1])
}
  print(hist(x,main=paste0("dist=",dist," n=",n," param=",
                           paste(as.character(param),collapse = ",")))) #put your code here
  }

hist_dist("exp",1000,param=1)
hist_dist("norm",1000,param=c(0,1))
hist_dist("bino",1000,param=c(100,0.5))
#ecc.


# 2 Calcolo della probabilità data una certa distribuzione####

## 2.1. Distribuzione esponenziale ####
load("data/campioni casuali/esponenziale.Rdata")
#standard plot di R
hist(x, freq = F, breaks = 20)
#generare la distribuzione teorica
x_g <- seq(0, 5, .1)
y <- exp(-x_g) #distribuzione di probabilità teorica
lines(x_g, y, col="red") # aggiungiamo la linea
dexp(1) # <- qual'è la probabilità di ottenere 1 da tale distribuzione?
abline(h=dexp(1),col="orange") #aggiungiamo oggetti grafici
abline(v=1,col="orange")
points(1,dexp(1),col="orange")

#funzione per plottare la distribuzione e il valore target
plot_prob_exp <- function(n,lambda,target){
  x <- rexp(n,lambda) #generiamo un campione di 100 elementi
  p <- dexp(target,lambda)
  hist(x, freq = F, breaks = n/5,
       main=paste0("n = ",n,", lambda = ",lambda,", valore target = ",target))
  x_g <- seq(0, 10, .1)
  y <- lambda*exp(-(lambda)*x_g) #distribuzione di probabilità teorica
  lines(x_g, y,col="red") # aggiungiamo la linea
  abline(h=p,col="orange") #1?
  abline(v=target,col="orange")
  points(target,p,col="orange")
  text(x=quantile(x,prob=.99),y=p+p/5,label=paste("p =",round(p,2)))
}
plot_prob_exp(100,1,2)
plot_prob_exp(100,5,.5)
plot_prob_exp(1000,.5,5)
plot_prob_exp(1000,.5,.5)

#using ggplot
ggplot_prob_exp <- function(n,lambda,target){
  x <- rexp(n,lambda) #generiamo un campione di 100 elementi
  x_g <- seq(0, 10, .1)
  y <- lambda*exp(-(lambda)*x_g) #distribuzione di probabilità teorica
  p <- dexp(target,lambda)
  ggp <- ggplot() +
    geom_histogram(aes(x = x, y = after_stat(density)),
                   bins = n/5,
                   col = "black") +
    geom_line(aes(y = y, x = x_g), col = "red")+
    geom_hline(yintercept = dexp(target,lambda),linetype=2,col="orange")+
    geom_vline(xintercept = target,linetype=2,col="orange")+
    geom_point(aes(x=target,y=dexp(target,lambda)),col="orange")+
    geom_text(aes(x=quantile(x,prob=.99),y=p+p/5,label=paste("p =",round(p,2))))
  return(ggp)
}

ggplot_prob_exp(100,1,2)
ggplot_prob_exp(100,5,.5)
ggplot_prob_exp(1000,.5,5)
ggplot_prob_exp(1000,.5,.5)

# la funzione dexp, pexp, qexp e rexp accettano più argomenti
# ad esempio la probabilità di ottenere 1 da tre distribuzioni esponenziali
# con
dexp(c(1,2),1)
lambda <- c(0.5, 1, 3)
x_g <- seq(0, 10, .1)
y<-c()
for (l in lambda) {
  y <- c(y,l*exp(-(l)*x_g))
}
df <- data.frame(y=y,
                 lambda=rep(lambda,each=101),
                 x=rep(x_g,3))

ggp <- ggplot(df)+
  geom_line(aes(y=y,x=x,col=as.factor(lambda)))+
  scale_color_discrete(name="lambda")
ggp

df$prob_1 <- dexp(1,df$lambda)
ggp + geom_vline(xintercept=1,linetype=2)+
  geom_point(data=df,aes(x=1,y=prob_1,col=as.factor(lambda)))

p_1 <- dexp(1,lambda) # qual'è il lambda a cui corrisponde la probabilità
print(paste("con lambda =",lambda," p(x=1) =", p_1))
# più alta?

# E con x=0.5?
df$prob_0.5 <- dexp(.5,df$lambda)
ggp + geom_vline(xintercept=.5,linetype=2)+
  geom_point(data=df,aes(x=.5,y=prob_0.5,col=as.factor(lambda)))

p_0.5 <- dexp(0.5,lambda) # qual'è il lambda a cui corrisponde la probabilità
print(paste("con lambda =",lambda," p(x=1) =", p_0.5))

# quantili
ggplot_prob_exp(1000,1,1)
print(paste("probabilità che il valore estratto sia uguale a 1 =",
            round(dexp(1),2))) # 
pexp(1) # probabilità che il valore sia minore o uguale di 1
pexp(1, lower.tail = F) # probabilità che il valore sia maggiore di 1
qexp(0.6321206)


## 2.2. Distribuzione Gaussiana (per casa) ####

## 2.3. Distribuzione di Bernoulli (per casa) ####


# 2. Il dataset AgrImOnIA ####
rm(list=ls())
dev.off()
gc()
# record page
# https://zenodo.org/records/7956006

# download
download.file("https://zenodo.org/records/7956006/files/Agrimonia_Dataset_v_3_0_0.Rdata?download=1",
              destfile = "data/Agrimonia_dataset.Rdata")

download.file("https://zenodo.org/records/7956006/files/Metadata_Agrimonia_v_3_0_0.csv?download=1",
              destfile = "data/Agrimonia_metadata.csv")

load("data/Agrimonia_dataset.Rdata")

Agrimonia_df <- AgrImOnIA_Dataset_v_3_0_0
rm(AgrImOnIA_Dataset_v_3_0_0)

names(Agrimonia_df)

metadata <- read.csv("data/Agrimonia_metadata.csv",sep = ";")
View(metadata)

head(Agrimonia_df)

length(unique(Agrimonia_df$IDStations))

subset_df <- subset(Agrimonia_df,IDStations=="1266")
summary(subset_df)
y<-subset_df$AQ_no2
y<-y[!is.na(y)]
mean(y)
quantile(y)
var(y)
hist(y,freq = F)
lines(density(y))

y<-log(y)
mean(exp(y))
quantile(exp(y))
var(exp(y))
hist(y,freq = F)
lines(density(y))

media <- mean(y)
exp(media)

x_g <- seq(min(y),max(y),.1)
y_teo <- (1/sqrt(2*pi*sd^2))*exp(-.5*((x_g-media)/sd)^2)
# abline(x_g,y_teo) ?

library(ggplot2)
ggplot()+
  geom_histogram(aes(x=y,y=after_stat(density)),col="black")+
  geom_line(aes(x=x_g,y=y_teo),col="orange",linewidth=2)

#ottenere la probabilità in relazione al valore
pnorm(log(10),media,sd)
pnorm(log(15),media,sd)
pnorm(log(20),media,sd)
pnorm(log(50),media,sd)
pnorm(log(50),media,sd,lower.tail = F)

#ottenere il valore in relazione alla probabilità
exp(qnorm(0.1,media,sd))
exp(qnorm(0.2,media,sd))
exp(qnorm(0.5,media,sd))
exp(qnorm(0.2,media,sd,lower.tail = F))



# install.packages("fitdistrplus")