#------------------------------------------------- # 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 #