# Primo esercizio

set.seed(2023)
albumina <- signif(sort(rnorm(65,mean=48,sd=5)),2)
albumina[1]<-albumina[1]-12
albu.b<-seq(24.5,65,by=5)

albumina

h.1<-hist(albumina,breaks=albu.b,xaxt="n",xlim = c(24,66),
     main="Albumina",xlab="[g/l]",ylab = "Frequenza")

axis(side=1,at=seq(0,66,by=1))
abline(v=seq(25,65),lty=3,col="grey")
abline(h=seq(1,30),lty=3,col="grey")

library(agricolae)

(t.1<-table.freq(h.1))
(m.pond <- sum(t.1$Main * t.1$Frequency)/sum(t.1$Frequency)) 

summary(albumina)

boxplot(albumina,horizontal = T)





# Importare il datasel

gamma <- read.csv("~/Scuola/SSMT/TAB-TRM/R2022/TAB2/papsat2223/gamma.csv", sep="")
str(gamma)
summary(gamma)
head(gamma)
gamma$gamma

# Istotogramma appropriato

hist(gamma$gamma)

# Dati arrotondati al primo decimale

ampl.gamma <- (max(gamma$gamma)-min(gamma$gamma))/(3+3.3*log10(length(gamma$gamma)))
ampl.gamma

ampl.gamma <- (max(gamma$gamma)-min(gamma$gamma))/(nclass.Sturges(gamma$gamma))
ampl.gamma


# Ampiezza 1, confine delle classi reali a x.x5 
gamma.b <- seq(min(gamma$gamma)-0.05,max(gamma$gamma)+0.05,by=1)
gamma.b

hist(gamma$gamma,breaks = gamma.b)

# Mediana e quantili, 95 percentile
quantile(gamma$gamma,c(0.025,0.05,0.25,0.5,0.75,0.95,0.975))

# Media
summary(gamma$gamma)

# Soglia dx 5% con ipotesi normale
qnorm(0.95,mean=mean(gamma$gamma),sd=sd(gamma$gamma))
boxplot(gamma$gamma,horizontal=T)

library(car)
qqPlot(gamma$gamma)

## Bonus
hist(gamma$gamma,breaks = gamma.b,probability=T)
gamma.d<-density(gamma$gamma)
lines(gamma.d,col="red")

library(spatstat)
quantile(gamma.d,0.95)
quantile(gamma$gamma,0.95)
qnorm(0.95,mean=mean(gamma$gamma),sd=sd(gamma$gamma))


# Regressione

regressione22 <- read.csv("~/Scuola/SSMT/TAB-TRM/R2022/TAB2/papsat2223/regressione22.csv")
regressione22 <- read.csv("regressione22.csv")

str(regressione22)

plot(regressione22$x,regressione22$y,xlim = c(0,10),ylim = c(0,15) )

fit <- lm(y ~ x, data = regressione22 )
abline(fit)
text(2,13,paste0("y = ",signif( coef(fit)[2],3),"*x + ",signif( coef(fit)[1])))
text(2,12,paste0("r² = ",signif(summary(fit)$r.squared,3)))

par(mfrow=c(1,2))
plot(fit,which=1:2)
par(mfrow=c(1,1))

## Bonus
#Non linerità eteroschedastica. Applicare logaritmo?

fit.l <- lm(log(y) ~ x, data = regressione22)
plot(regressione22$x,log(regressione22$y))

par(mfrow=c(1,2))
plot(fit.l,which=1:2)
par(mfrow=c(1,1))

abline(fit.l)
text(5,2.6,paste0("y = ",signif( coef(fit.l)[2],3),"*x + ",signif( coef(fit.l)[1])))
text(5,2.4,paste0("r² = ",signif(summary(fit.l)$r.squared,3)))

par(mfrow=c(1,2))
plot(fit,which=1:2)
par(mfrow=c(1,1))