#########################################################
# Parte 1. : calcolo  intervallo  di confidenza
#            utile per fare esperimenti varian-
#            do il valore di  parti delle equa-
#            zioni ed indagare l'effetto sullo
#            intervallo di confidenza ottenuto.
#########################################################

#Campione casuale di n=11 donne scelte da una popolazione omogenea dove si e' misurata la quantita' giornaliera di energia assunta (in kjoule).

energia <- c(5260, 5470, 5640, 6180, 6390, 6515, 6805, 7515, 7515, 8230, 8770)

# verifica visuale normalita' 
qqnorm(energia)
qqline(energia)

# test di normalita'
shapiro.test(energia)


# Assumendo probabilita' normale costruiamo l'intervallo di confidenza per il quantitativo medio di energia assunta giornalmente nella popolazione.
# Supponiamo di *SAPERE* che la deviazione standard della popolazione sia sigma=1000 kj . Allora l'intervallo di confidenza al 95% per mu e':

#  xbar  +- 1.96sigma / sqrt(n)

# ossia un intorno della media campionaria

mean(energia)

# di raggio (il cosiddetto margine di errore) 1.96 * 1000/sqrt(11) . QUindi l'intervallo di confidenza al 95% per mu ha estremi:

mean(energia) - 1.96 * 1000/sqrt(length(energia))
mean(energia) + 1.96 * 1000/sqrt(length(energia))


####### CASO IN CUI SIGMA NON E' NOTO #########

# In questo caso l'intervallo di confidenza al 95% e' : 

#  xbar +- t * s/sqrt(n)
# Dato che s (la deviazione standard campionaria) / sqrt(n) e' l'errore standard (SE) posso scrivere equivalentemente che in questo caso l'intervallo di confidenza e' :
 #  xbar +- t * SE
 xbar  -  sigma*1.96/sqrt(length(energia))

# dove:
# s e' la deviazione standard del campione (non quella supposta della popolazione)
# t e' un coefficiente maggiore di 1.96 che dipende dalla dimensione del campione

# Il coefficiente t si ottiene non dalla tavola della normale ma da quella della distribuzione t di Student. Si puo' sostituire la tavola della t con la seguente funzione R (qt):

LC <- 0.95  # il livello di confidenza richiesto
n <- 11     # numerosita' del campione
t <- qt( LC + (1-LC)/2, n-1)
t

# qt() e' l'equivalente t di Student della qnorm() della normale. Anche per t di Student esistono dt()  (equivalente di dnorm) e pt (eqiovalente di pnorm).

# Ora che abbiamo calcolato t possiamo calcolare l'intervallo di confidenza al 95% anche nel caso in cui il valore di sigma della popolazione non e' dipsonibile. In R si fa cosi' :

SE <- sd(energia)/sqrt(length(energia))
xbar <- mean(energia)
# calcolo estremi intervallo:
A <- xbar - t * SE
B <- xbar + t * SE
c(A,B)

# Notate che e' piu' ampio del precedente...
# Funzione per il calcolo diretto dell'intervallo di confidenza usando il metodo basato su t (sigma popolazione non nota):

t.test(energia)
t.test(energia, conf.level=0.99)

######################################################################
# Parte 2. : simulazione di campioni per indagine su C.I. 
#
# Quanto "casuali" sono i CI ? Cosa e' il livello di "confidenza"? 
#
#          utile per capire come e' possibile simulare
#          campionamenti per  studiare la procedura di
#          calcolo di C.I. e p-value
#####################################################################
#          Obiettivi:
#                Informatica:     - procedura di simulazione
#                                 - utilizzo di user defined function
#                                 - conoscenza funzione apply()
#                                 - conoscenza strutture dati
#                Statistica:
#                                 - interpretazione dei risultati della
simulazione
#####################################################################

# Dataset: lunghezza coda lucertole (n=24) ##########################

lizard <- c(6.2, 6.6, 7.1, 7.4, 7.6, 7.9, 8, 8.3, 8.4, 8.5, 8.6, 8.8, 8.8, 9.1, 9.2, 9.4, 9.4, 9.7, 9.9, 10.2, 10.4, 10.8, 11.3, 11.9)


# Imposto numero di campionamenti (dalla normale) ##################
n.draw <- 100

# Voglio campionare da una normale con una mu *DECISA* a priori ( poi vedremo se i C.I. ottenuti contengono il valore di mu utilizzato.
mu <- 9
# Dimensione dei campioni estratti (campionati dalla normale via rnorm()
n <- 24

# Come SD usiamo la SD campionaria (ottenuta dal campione originale, il
vettore lizard) 
SD <- sd(lizard)

# STEP I: creo matrice contenente 100 vettori (campioni simulati estratti
dalla normale via rnorm) ognuno di numerosita' n=24 (come il vettore lizard
contenente i dati originali). I 100 vettori costituiscono le COLONNE della
matrice che, quindi, ha dimensione 24 x 100. Salvo matrice in variabile draws.

draws <- matrix(rnorm(n.draw * n, mu, SD), n)

# STEP II: creo funzione per estrarre GLI ESTREMI dell'intervallo di
confidenza ottenuto mediante un t.test(). La funzione, dato un vettore di 24
osservazioni (ossia UNA COLONNA della matrice draws), ritorna un vettore di 2
elementi contenente il limite INFERIORE ed il limite SUPERIORE dell'intervallo
di confidenza calcolato per la media della popolazione. Chiamo la funzione
get.conf.int() 

get.conf.int <- function(x){ t.test(x)$conf.int }

# STEP III: calcolo intervallo di confidenza per tutti i campioni simulati
(estratti dalla normale) ossia per TUTTE le colonne della matrice draws.
Quindi otterro' un numero di C.I. pari al numero di colonne della matrice
della simulazione (quindi 100). Per farlo utilizzo funzione apply() applicando
la funzione creata in step II a tutte le colonne (dimensione 2) della matrice
draws. Salvo tutto nella matrice conf.int

conf.int <- apply(draws, 2, get.conf.int)

# STEP IV: conto quante volte (sulle 100 simulazioni effettuate) la vera mu
della popolazione (che conosco perche' l'ho scelta io mentre creavo i
campionamenti, ed era stata impostata a 9) cade nel C.I. calcolato. 

sum(conf.int[1, ] <= mu & conf.int[2, ] >= mu)

# STEP V: creo grafico con tutti i C.I. calcolati . Le fasi sono 3. a) creo
grafico vuoto. b) disegno tutti i C.I. c) disegno riga linea vertivale in
corrispondenza del valore di mu utilizzato in fase di simulazione dei campioni
(ossia mu=9). E' possibile notare che non tutti i C.I. calcolati contengono
mu. 


plot(range(conf.int), c(0, 1 + n.draw), type = "n", xlab = "mean tail length", ylab = "sample run")
for (i in 1:n.draw){ lines(conf.int[, i], rep(i, 2), lwd = 2)  }
abline(v = 9, lwd = 2, lty = 2)


############## ESERCIZIO #####################################################
# Nel test con le simulazioni il livello di confidenza utilizzato dalla
# funzione e' quello di default, 0.95. Modificate la funzione per passare
# voi il valore desiderato. (leggete il manuale di t.test() e fate attenzione
# al parametro conf.level