#-------------------------------------------------
# analisi dati sulla misura delle densita'
# di oggetti irregolare ('caramelline di vetro'),
# con spunti di riflessione sui 'fit'
#
#     G. D'Agostini, nov 09
#-----------------------------------------------------

# 0)  ---------------------------------------------------
# leggiamo, guardiamo e plottiamo i dati
dati <- read.table("dati_densita.txt")
dati
volume <- dati$V1
massa  <- dati$V2
n      <- length(massa)

# prima stima della densita':
rho.0 <- (massa[n]-massa[1])/(volume[n]-volume[1])


plot(volume, massa)

# 1)  ------------------------------------------------- 
# fit massa Vs volume
fit.mv <- lm(massa ~ volume)
fit.mv


# sovrapponiamo il fit sul grafico
abline(fit.mv, col='blue')


# il risultato di lm() e' una lista con tante informazioni
str(fit.mv)

# il cui 'summary' e' dato da
summary(fit.mv)

# per ora siamo interessati a
fit.mv$coefficients
# ( si posono visualizzare anche con
coef(fit.mv)
# )

# per ripulirli dai nomi:
as.vector(fit.mv$coefficients)

# il secondo dei coefficienti ci da' la densita',
# come coefficiente di linearita' fra massa e volume
rho <- as.vector(fit.mv$coefficients)[2]

# 2)  -----------------------------------------------------
# Confrontiamo lm() con il metodo dei 'minimi quadrati'
#  
# funzione che calcola la somma dei quadrati degli scarti
sq <- function(x) {
  # x[1]: intercetta
  # x[2]: 'pendenza' (coeff. di linearita')
  scarti <- massa - (x[1] + x[2] * volume)
  return( sum( scarti^2) )
}

fit.mv.mq <- nlm(sq, c(0, 1),steptol=1e-10, gradtol=1e-10 )
# attenzione nlm e' molto sensibile al punto iniziale!!
# con steptol e gradtol di default non funziona bene!!
#--------------------------------------------------------------------------


# 3) -----------------------------------------------------
# somma moduli scarti
# funzione che calcola la somma dei moduli degli scarti
sms <- function(x) {
  # x[1]: intercetta
  # x[2]: 'pendenza' (coeff. di linearita')
  #cat(sprintf("c, m: %g,  %g \n", x[1], x[2]))
  scarti <- massa - (x[1] + x[2] * volume)
  #print(scarti)
  s <- sum( abs(scarti) )
  return(s)
  #print(s)
#  if (is.finite(s)) {
#    return(s)
#  } else {
#    #
#    cat(sprintf("ATT: c, m: %g,  %g \n", x[1], x[2]))
#    print(scarti)
#    return(1e+100)
  }
}

# proviamo questa, partendo dal risultato precedente
# (att. nlm sembra alquanto instabile...)
fit.mv.sms <- nlm(sms, c(fit.mv.mq$estimate[1], fit.mv.mq$estimate[2]),
                  steptol=1e-10, gradtol=1e-10 )

# risultati
fit.mv.sms$estimate

# sembrano 'molto simili': 
# (ma, in realta' per decidere veramente
# quanto sono simili bisogna imparare a capire
# anche con quante cifre prendere il risultato sulla
# densita'

#4) minimizziamo somma delle distanze retta-punti
#  -> lasciato come esercizio (se si pensa che possa avere un senso)
#  -> si raccomanda di far partire nlm() dalla soluzione
#     dei minimi quadrati, come fatto sopra


# 5) metodo 'matriciale' ---------------------------------------------
# costruiamo questa 'strana' matrice ("design matrix"):
X <- matrix( rep(1,n), n, 1)
X <- cbind(X, volume)
X

t.X       <- t(X)                       # trasposta
t.X.X     <- t.X %*% X                  # prodotto fra t(X) e X
i.t.X.X   <- solve(t.X.X)               # inversione 
parametri <- i.t.X.X %*% t.X %*% massa  # <= risultato
parametri

# in forma piu' compatta,
solve(t(X) %*% X) %*% t(X) %*% massa  

# confrontiamo con risultato di lm():
(parametri[2]-rho)/rho

# sono praticamente la stessa cosa!
# -> quindi probabilmente lm() internamente usa questo metodo


#------------------------------------------------------------
# 6) analisi dei residui
# analiziamo i 'residui'
fit.mv$residuals
plot(dati$V1, fit.mv$residuals)
abline(0,0, col='red')

# confrontiamo con
residui <- massa - X %*% parametri
residui - fit.mv$residuals
# effettivamente sono la stessa cosa!

# istogramma dei residui:
hist(residui)
# ... e loro media
mean(residui)  # praticamente zero


# interessanti a questo punto
# - la media dei moduli
mean(abs(residui))
# - la deviazione standard
sd(residui)
# -> che significato hanno?


# ... guardiamo ad altre informazioni fornite da lm()
summary(fit.mv)
# -> vengono forniti anche errori sui coefficienti
# -> occhio al "Residual standard error",
#    prossimo a sd(residui)

# proviamo a calcolare:
sd (residui) * sqrt((n-1)/(n-2))
# e
sqrt(sum(residui^2)/(n-2))
# -> quindi il "Residual standard error" e'
#    esattamente questa 'cosa' qui. Perche?
# -> Nota: che dimensioni fisiche e unita' di
#          misura hanno queste quantita'? 

#---------------------------------------------------------
# 7) fit lm() scambiando le variabili:

plot(massa, volume)

# 1)  ------------------------------------------------- 
# fit massa Vs volume
fit.vm <- lm(volume ~ massa)
fit.vm


# sovrapponiamo il fit sul grafico
abline(fit.vm, col='blue')

#chiaramente
rho.vm <- 1/as.vector(fit.vm$coefficients)[2]
# confrontiamo:
(rho.vm-rho)/rho
# molto molto simili, 
# MA non proprio identici entro le approssimazioni numeriche!

# cosa significano in questo caso 
fit.vm$residuals
# e 
sqrt(sum(fit.vm$residuals^2)/(n-2)) # ?
# Che dimensioni hanno?
# Si noti come quest'ultimo e' uguale, nuovamente,
# al "Residual standard error", fornito da
summary(fit.vm)

#-------------------------------------------------------------
# 8) altri studi suggeriti
#
# - aumentare artificialmente gli errori, ad esempio con
#       volume.sim <- volume + sv*rnorm(n)
#   oppure
#       massa.sim <- massa + sm*rnorm(n)
#   ove sv e sm sono numeri a piacere, inizialmente 1,
#   che danno la scala degli errori
#
# - provare a confrontare i risultati dei fit Massa Vs volome
#   e volume Vs massa a mano a mano che gli errori crescono
#