library(ggplot2)
set.seed(1)
# dallo script Statistica_3alezione: campionaria vs teorica ####
#funzione esponenziale
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)
  }
}

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, .01))) {
    x <- rexp(n, lambda)
    df <- rbind(
      df,
      data.frame(
        n = n,
        lambda = lambda,
        media_teo = 1 / lambda,
        media_camp = mean(x),
        var_teo = 1 / lambda ^ 2,
        var_camp = var(x)
      )
    )
  }
}

View(df)

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

p <-
  p + scale_x_continuous(breaks = 10 ^ c(1, 2, 3, 4, 5, 6), trans = "pseudo_log")

p #notare l'asse delle X !

p + scale_y_continuous(trans = "reciprocal") # y_new = 1/y_old invertito l'ordine!

# DOMANDA: come si distribuisce la media?
# RISPOSTA: 
# suggerimento: vedi il Teorema del limite centrale

hist(df$media_camp[df$lambda == .5], freq = F, ylim = c(0, 3))
lines(density(df$media_camp[df$lambda == .5], width = .5)) # check width parameter

media_delle_medie <- mean(df$media_camp[df$lambda == .5])
sd_delle_medie <- sqrt(var(df$media_camp[df$lambda == .5]))

abline(v = media_delle_medie, col = "red")
abline(v = media_delle_medie - sd_delle_medie, col = "yellow")
abline(v = media_delle_medie + sd_delle_medie, col = "yellow")
abline(v = media_delle_medie - 1.96 * sd_delle_medie, col = "orange")
abline(v = media_delle_medie + 1.96 * sd_delle_medie, col = "orange")

text(
  y = 2.5,
  x = media_delle_medie - .05,
  labels = expression(hat(mu)),
  col = "red"
)

text(
  y = 2.5,
  x = media_delle_medie - sd_delle_medie - .05,
  labels = expression(hat(sigma)),
  col = "yellow"
)
text(
  y = 2.5,
  x = media_delle_medie + sd_delle_medie + .05,
  labels = expression(hat(sigma)),
  col = "yellow"
)

text(
  y = 2.5,
  x = media_delle_medie - 1.96 * sd_delle_medie - .09,
  labels = "2.5th perc",
  col = "orange",
  srt = 90
)
text(
  y = 2.5,
  x = media_delle_medie + 1.96 * sd_delle_medie + .09,
  labels = "97.5th perc",
  col = "orange",
  srt = 90
)

#vediamo se gli intervalli di confidenza teorici assomigliano a quelli empirici
perc_2.5 <- media_delle_medie - 1.96 * sd_delle_medie
perc_97.5 <- media_delle_medie + 1.96 * sd_delle_medie

medie <- df$media_camp[df$lambda == .5]
medie <- medie[order(medie)]
medie
0.025 * length(medie)
0.975 * length(medie)
medie[0.025 * length(medie)] #quanto c'è di differenza con il percentile teorico?
# risposta:
medie[0.975 * length(medie)] #quanto c'è di differenza con il percentile teorico?
# risposta:

# DOMANDA: Qual'è la varianza teorica della media?
# cioè VAR[1/n*sum(x)] = VAR[x]/n
# RISPOSTA: VAR[x]/n

# per ogni media campionaria possiamo calcolare la sua varianza
df$var_N <- df$var_camp / df$n
df$sd_N <- sqrt(df$var_N)

p + geom_ribbon(
  aes(
    ymin = df$media_camp - 1.96 * df$sd_N,
    ymax = df$media_camp + 1.96 * df$sd_N,
    x = n,
    col = as.factor(lambda),
    fill = as.factor(lambda),
  ),
  alpha = 0.4,
  linetype = 3
) + theme_light()

# verifichiamo che la media teorica cada dentro gli intervalli di confidenza
# della media campionaria il 95 % circa delle volte

df$CI_lwr <-
  df$media_camp - 1.96 * df$sd_N #definiamo il limite inferiore
df$CI_upp <-
  df$media_camp + 1.96 * df$sd_N #definiamo il limite superiore
sub_df <- subset(df, lambda == 0.5) #guardiamo solo ad un lambda

#quante volte la media teorica è fuori dagli intervalli di confidenza?
sum(sub_df$media_teo > sub_df$CI_lwr &
      sub_df$media_teo < sub_df$CI_upp) / nrow(sub_df)
# DOMANDA: di quanto cambia rispetto al valore teorico?
# RISPOSTA:


# Tiro dei dadi ####
set.seed(1)
rm(list = ls())
gc()
n <- 100000 #numero di tiri
facce <- 10 #numero di facce
tiri <- sample(1:facce, n, replace = T)
plot(table(tiri))
df <- as.data.frame(table(tiri))
df <- merge(data.frame(tiri = 1:6), df, all = T)
df$Freq[is.na(df$Freq)] <- 0
ggplot(df) +
  geom_col(aes(y = Freq, x = 1:facce))

# DOMANDA: come si calcolano le frequenze relative?
# RISPOSTA:

df$Freq_rel <- df$Freq / n
p <- ggplot(df) +
  geom_col(aes(y = Freq_rel, x = tiri))
p

# Stiamo considerando una distribuzione discreta uniforme
# DOMANDA: come si calcolano le frequenze relative teoriche?
# RISPOSTA:

freq_teo <- 1 / facce
p <- p + geom_col(aes(y = freq_teo, x = 1:facce),
                  col = "red",
                  fill = NA)
p

#aggiungiamo le etichette
p <- p + geom_text(aes(
  x = 3,
  y = max(Freq_rel),
  label = "RED=theoretical
BLACK=empirical"
))
p

#simuliamo 10000 campioni con n crescenti
rm(list = ls())
gc()
facce <- 6
media_camp <- rep(NA, length(seq(10, 10 ^ 5, 10)))
var_camp <- rep(NA, length(seq(10, 10 ^ 5, 10)))
for (n in seq(10, 10 ^ 5, 10)) {
  tiri <- sample(1:facce, n, replace = T)
  idx <- which(seq(10, 10 ^ 5, 10) == n)
  media_camp[idx] <- mean(tiri)
  var_camp[idx] <- var(tiri)
}

p_mean <- ggplot() +
  geom_line(aes(y = media_camp, x = 1:length(media_camp))) +
  scale_x_continuous(breaks = c(10 ^ c(1, 2, 3, 4, 5)), trans = "pseudo_log") +
  ggtitle(label = "sample mean for increasing number of",
          subtitle = "increasing number of dice rolls") +
  xlab("numero di lanci") +
  ylab("media campionaria")
p_mean

p_var <- ggplot() +
  geom_line(aes(y = var_camp, x = 1:length(media_camp))) +
  scale_x_continuous(breaks = c(10 ^ c(1, 2, 3, 4, 5)), trans = "pseudo_log") +
  ggtitle(label = "sample variance for increasing number",
          subtitle = "of increasing number of dice rolls") +
  xlab("numero di lanci") +
  ylab("varianza campionaria")
p_var

# DOMANDA: qual'è la media di una distribuzione uniforme discreta?
# RISPOSTA:

# suggerimento: formula somma integrali -> sum(1:n) = n(n+1)/2
# suggerimento: E[1/n*sum(1:n)] = (n+1)/2

media_teo <- (facce + 1) / 2

p_mean <- p_mean + geom_hline(yintercept = media_teo,
                              col = "red",
                              linetype = 2)
p_mean

# DOMANDA: qual'è la varianza di una distribuzione uniforme discreta?
# RISPOSTA:
# suggerimento: formula somma quadrati -> sum(x^2) = n(n+1)(2n+1)/6
# suggerimento: E[(x - E[x])^2] = E[x^2] - E[x]^2 = 1/n*sum(x^2) - ((n+1)/2)^2

var_teo <- (facce + 1) * (facce - 1) / 12
p_var <- p_var + geom_hline(yintercept = var_teo,
                            col = "red",
                            linetype = 2)
p_var

# generiamo m campioni di n dati da un uniforme discreta
# (tiro dei dadi)
rm(list = ls())
gc()
facce <- 6
campione <- rep(10 ^ c(1, 2, 3, 4, 5),each = 1000)
df <- data.frame(matrix(NA,ncol = 3,
                        nrow = length(campione)))
names(df)<-c("media_campionaria","varianza_campionaria","n_campione")
i<-0
for (n in campione) {
  tiri <- sample(1:facce, n, replace = T)
  i <- i + 1
  df$media_campionaria[i] <- mean(tiri)
  df$varianza_campionaria[i] <- var(tiri)
  df$n_campione[i] <- n
}

df$n_campione <- as.factor(df$n_campione)

summary(df)
df_10 <- subset(df,n_campione==10)
media_10 <- df_10$media_campionaria
hist(media_10)
hist(media_10,freq = F, ylim=c(0,0.8))
lines(density(media_10))
lines(density(media_10,width = 1),col="red")
lines(density(media_10,width = 2),col="orange")
lines(density(media_10,width = 3),col="yellow")
lines(density(rnorm(1000,mean(media_10),
                    sqrt(var(media_10)))),col="purple")
for (i in 1:1000) {
  lines(density(rnorm(1000,mean(media_10),
                      sqrt(var(media_10)))),col="purple")
}

media_medie <- mean(media_10)
sd_medie <- sqrt(var(media_10))
media_medie
sd_medie
normale_ci <- 1.96*sd_medie
upp <- media_medie + normale_ci
lwr <- media_medie - normale_ci
sum(media_10 < lwr)
sum(media_10 > upp)

#run it over time!

#Il test t di Student per la verifica d'ipotesi su un valore medio
#t.test()
mean(media_10)
t.test(media_10)
t.test(media_10,alternative = "greater")
t.test(media_10,alternative = "less")
t.test(media_10,mu = 2)
t.test(media_10,mu = 2,alternative = "greater")
t.test(media_10,mu = 2,alternative = "two.sided", conf.level = .1)
t.test(media_10, mu = mean(media_10),conf.level = .95)
t.test(media_10, conf.level = .95)
hist(media_10)
# sono diversi dai quantili della distribuzione
quantile(media_10, probs = c(.025,.975))
media_10[order(media_10)][.025*length(media_10)]
media_10[order(media_10)][.975*length(media_10)]

# DOMANDA: come standardizzare i valori?
# RISPOSTA:
media_10_std <- (media_10 - mean(media_10)) / sqrt(var(media_10))
hist(media_10_std)
t.test(media_10_std)