##esercizio 0
#0 se non piove, 1 se piove per la var Xi,
#prove indipendenti
#p è la probabilità che piova o meno ed è tra 0 e 1

#1. Scrivere in funzione di p la plegge di probabilità di Xi per un generico i
#Dato che chiede di UN generico i, do la definiizone di legge di probabilità della singola osservazione
#quindi dato che ogni osservazione è una bernoulliana, fx(i)=p(1-p)

#2. indicando rispettivamente E1, E2 gli eventi "piove in almeno un giorno infrasettimanale"
#e "non piove nel fine settimana"
#foglio

#3. gli eventi E1, e2 sono indipendenti? Sono anche mutualmente esclusivi?
#si sono indipendenti, perchè gli eventi sono indipendenti tra loro data l'ipotersi che
#le condizioni metereologiche di un giorno non influenziano le condizioni meteo dei giorni rimanenti
#Non possono considerare la mutua esclusione, dato che questa prevede che il verificarsi di uno, impedisca il verificarsi dell'altro
#Quindi se E1 ha una probabilità p, E2 deve avere una probabilità==0

#4. sapendo p=0.4, calcola la probabilità di
#E1, E2, E3=piove almeno in un giorno infrasettimanale, ma non durante il finesettiman"
#E4="piove almeno in un giorno infrasettimanale OPPURE non piove nel fine settimana"
p=0.4
E1=1-(pbinom(0,5,p));E1
E2=1-(pbinom(0,2,(1-p))); E2
pbinom(2,2,p)
E3=E1*E2; E3
E4=E1+E2-E3; E4

##esercizio 1
#Indichiamo B la var aleatoria che conta il numero di giornate piovese in una settimana
#1. Punti di massa di B?
#i punti di massa sono i valori da 0 a 7 compresi gli estremi

#2. scruvere in funzione di p la legge di probabilità di B
#uso la funzione di massa di probabilità della binomiale foglio

#3. sono sufficienti 7 giorni di osservazioni per poter affermare che è maggiore di 0.4 la probabilità
#di compiere nella stima p un errore invferiore a 0.25 (si ricordi p(1-p) non può mai superare il valore di 1/4)
#usare chebyshev perchè n è finito e non tendente a infinito

##esercizio 2
#inserire il dataframe
dir()
dd=read.delim(file="simulazione.csv", header=TRUE, sep=";")
attach(dd)
str(dd)
#1. quante centraline ci sono in lombardia:
#considerando che ogni riga è una centralina e tutte le province sono lombarde, mi basta contare
#il numero di righe
numeroCentraline=nrow(dd); numeroCentraline
m=nrow(dd)
#2. come sono distribuite le centraline nelle diverse province della lombardia
barplot(table(PROVINCIA))
#si nota il picco di centraline (circa 35) nella provincia di milano

#3. Quale provincia è meno rappresentata? Qual è la pià rappresentata? Como è la meno rappresentaa (meno di 5 centraline)
#invece milano è la più rappresentata (35 centraline)

#4. calcolare l'indice di eterogeneità per il carattere PROVINCIA e commentare il risultato ottenuto
gini=sum(prop.table(table(PROVINCIA))^2)*(m/(m-1))

#5. per ciascuno di essi, indicare in tabella1 medianda dev standard e distanza interquartile tra i 3 param
mean(C6H6, na.rm=TRUE)
mean(SO2,  na.rm=TRUE)
mean(CO,  na.rm=TRUE)
median(C6H6,  na.rm=TRUE)
median(SO2,  na.rm=TRUE)
median(CO,  na.rm=TRUE)
sd(C6H6,  na.rm=TRUE)
sd(SO2,  na.rm=TRUE)
sd(CO, na.rm=TRUE)
IQR(C6H6,  na.rm=TRUE)
IQR(SO2,  na.rm=TRUE)
IQR(CO,  na.rm=TRUE)
help(plot)

#6. mettere a confronto i 3 inquinanti con l'opportuna modalità grafica commentando il grafico ottenuto
x=c(0,8)
y=c(0,15)
plot(table(C6H6), type="b", xlim=x, ylim=y)
par(new=TRUE)
plot(table(SO2), type="b", xlim=x, col="green", ylim=y)
par(new=TRUE)
plot(table(CO), type="b", xlim=x, col="red", ylim=y)
#C6H6 -> si hanno poche rilevazioni e di queste, è stato registrato in più centraline un valore di 0.9
#il valore massimo raggiunto è 4, ma in una sola centralina
#SO2 -> in molte centraline è stato rilevato il valore 3, ha raggiunto il picco a 8 in una sola centralina
#CO -> in tantissime centraline(13) è stato registrato un valore di 0.5, ma il picco massimo è di 1.3 in una centralina

#Fagli elementi cadmio e bezoapirene abbiamo perso i dati, ma ci sono rimasti i grafici
#della funzione cumulativa empirica e della freqeuza relativa
#7. Completare la figura con le corrispondeze appropriate
#il primo ecdf è sicuramente il cadmio, perchè ho un grande salto iniziale, poi un piccolo salto
#un altro salto più piccolo del primo, ma sempre grande, infine salti di dimensione decrescente.
#
#nel secondo ecdf è sicuramente il benzoapirene, perchè si verifica un salto molto piccolo, seguito da uno grande,
#poi un'altro grande, ma minore del primo, uno piccolo, mancanza di osservazioni, poi altri 2 piccoli salti

#Consideriamo la concentrazione di ozono O3:
#8. calcolare la tabella dell frequnze relative della concentrazione annuale
str(dd)
prop.table(table(O3))

#9. Tracciare il grafico della funzione cumulativa empirica di O3
plot(prop.table(table(O3)))
plot(ecdf(O3))
#la ecdf rispetta il grfico delle frequenze empiriche

#10. quante centraline hanno registrato una concentrazione uguale a 46 mg/mcubo di O3?
d1=subset(dd,O3==46)
sum(d1[,c("O3")])/46 #4 centraline hanno rilevato 46 di O3, nrow restituisce null, nvm
nrow(d1) #4 uguale

#oppure
d1=O3[O3==46 & !is.na(O3)]
length(d1)

#11. quante centraline hanno registrato una concentrazione compresa tra 45 e 50 mg/mcubo compresi gli estremi di O3?
d2=subset(dd,O3>44 & O3<51)
nrow(d2) #25 osservazioni/centraline hanno rilevato una concnetrazione compresa tra 45 e 50 di O3

#12. Indicar quali valori di O3 sono "particolarmente elevati risptto agli altri, motivando con un grafico
boxplot(O3)
#per la provincia di lecco c'è una osservazione >90, precisamente 93, l'unico oulier presente
#si può ipotizzare un errore di battitura o un errore di inserimento. In ogni caso quella osservazione supera
#di gran lunga il range interquartile "superiore" del boxplot

#Consideriamo la variabile casuale X="numero di giorni in un anno nei quali viene superata la sogli di tolleranza di O3"
#abbiamo a disposizione un campione csuale X1,..,Xn di taglia n=150 estratto dalla popolazione X (la centralina i-esima
#è la variabile Xi) e memorizzato nel file scaricato
#13. si fornica una stima del valore atteso e della deviazione standard di X
#prima guardo se la soglia segue una qualche legge
str(dd)
barplot(O3_GIORNI_SUPERAMENTO_TOLLERANZA)
#è a campana, ma non sembra seguire una legge in particolare, stimo i due parametri indistintamente
#stimo il valore atteso tramite lo stimatore Tn=media campionaria

e=mean(O3_GIORNI_SUPERAMENTO_TOLLERANZA, na.rm=TRUE); e #il valore atteso di giorni in cui si supera la tolleranza è di 11,20
#stimo la varianza con lo stimatore Tn=varianza campionaria
var=var(O3_GIORNI_SUPERAMENTO_TOLLERANZA, na.rm=TRUE); var
help(var) #dato che per default su R la funzione var è la varianza campionaria, lo uso come stimatore
sd(O3_GIORNI_SUPERAMENTO_TOLLERANZA, na.rm=TRUE)#FUNZIONA perchè usa la varianza campionaria
sqrt(var)#della varianza stimata con la varianza campionaria, ci faccio la radice per ottenere la dev standard

#14. Ci sono buoni motivi per pensare che X segua la legge binomiale?  Si giustifichi la risposta
#e(x)=n*p
#var(x)=n*p(1-p)
#posso ricavare i p
#p=e(x)/n
p=e/150; p
x=c(0,30)
r=rbinom(150,150,p)
qqplot(O3_GIORNI_SUPERAMENTO_TOLLERANZA,r, pch=20, xlim=x, ylim=x)
par(new=TRUE)
qqplot(r,r,col="red", pch=20, xlim=x, ylim=x) #non so fare la qqline per la binom :cccc
#Non può seguire una legge binomiale, dato che la distribuzione è molto lontana dalla bisettrice