#-------------------------------------------------------
#
# analisi densita' delle "caramelline" di vetro
#
# GdA  PED, 17 dic 2010
#-------------------------------------------------------


m <- c(0, 6.6, 13.2, 33.6, 53.8)
v <- c(68.1, 71.0, 73.9, 82.0, 90.3)

fit <- lm(v~m)
str(fit)
plot(m,v)
abline(fit)
1/fit$coefficients[2]

#---------------------------------------------------
# conti ripetuti "a mano":
n <- length(m)
cov.mv <- mean(m*v) - mean(m)*mean(v) 
var.m  <- mean(m*m) - mean(m)^2
E.slope  <- cov.mv / var.m
E.interc <- mean(v) - slope * mean(m)
#--------------------------------------------------

sigma.v <- sqrt( sum(fit$residuals^2)/n )   # o .../(n-2) ->  provare
sigma.slope  <- sigma.v/sqrt(n) / sqrt(var.m)
sigma.interc <-  sqrt(mean(m^2))*sigma.slope
correlazione.slope.interc <- -mean(m)/sqrt(mean(m^2))

# incertezza relativa sulla pendenza
r.slope <- sigma.slope/E.slope

# da cui:
E.rho <- 1/E.slope
sigma.rho <- E.rho * r.slope

#-------------------------------------------------------------
# Cosa succede se si fa il fit massa in funzione del volume?
# [ovvero:  fit.inv <- lm(m ~ v) ] 
# In particolare, come cambia il risultato finale?
#
#-------------------------------------------------------------
# Provare a fare una function, funzione delle generiche x e y,
# fa il fit con lm() e calcola anche le incertezze, usando il metodo
# dei residui (eventualmente si puo' mettere nella chiamata
# un parametro per desidere se usare "n" o "n-2" 
#
# Questa potrebbe esserne l'ossatura:
lfit <- function(x,y) {
  ris <- NULL

  if (length(x) != length(x)) {
    ris$errore <- "Errore: x e y devono avere le stesse dimensioni"
    return(ris)
  }

  fit <- lm(y~x)

  ris$coefficients <- fit$coefficients
  # calcoliamo poi le sigma dei coefficienti
  # ....
  # ....
  # quindi
  ris$sigma.y <- sigma.y
  ris$sigmas  <- c(sigma.interc, sigma.slope)
  return(ris)
}