R/linearfit.R

In Fioroni1863179/LabMec: Prova Individuale LabMec

#' Fit lineare su una serie di misure ed estrapolazione di punti non misurati
#'
#' La funzione calcola la retta che meglio approssima i dati passati stampando coefficiente angolare, intercetta e coefficiente di correlazione. Se viene passato il vettore estX, per ogni valore del vettore calcola la migliore stima con incertezza.
#'
#' @param x
#' Il vettore contenente le ascisse dei punti misurati
#' @param y
#' Il vettore contenente le ordinate dei punti misurati
#' @param sigma
#' Il vettore contenente le incertezze sulle misure di \emph{y}. È accettato anche un singolo valore.
#' @param estX
#' DEFAULT = NULL
#'
#' Se il valore è diverso da NULL, per ogni elemento del vettore estrae l'ordinata con la rispettiva incertezza
#' @param verbose
#' DEFAULT = FALSE
#'
#' Se il valore è TRUE, durante i calcoli stampa sul terminale i valori di x_medio, y_mendio, xy_medio, x^2_medio, varianza_x, covarianza_xy, sommatoria di 1/sigma^2
linearfit <- function(x, y, sigma, verbose=FALSE, estX = NULL){
  if(missing(x)) stop("manca il vettore delle x")
  if(missing(y)) stop("manca il vettore delle y")
  if(missing(sigma)) stop("manca l'incertezza sulle y")

  if(length(sigma) == 1 ){
    sigma = rep(sigma, length(x))
  }
  invsigma <- 1/sigma^2
  meanx <- weighted.mean(x,invsigma)
  meany <- weighted.mean(y,invsigma)
  variance <- weighted.mean((x-meanx)^2,invsigma)
  Em <- (weighted.mean(x*y, invsigma) - meanx*meany)/variance
  sigm <- sqrt(1/(sum(invsigma)*variance))
  Ec <- meany - Em*meanx
  sigc <- sqrt(weighted.mean(x^2,invsigma)/(sum(invsigma)*variance))
  correl <- -meanx/sqrt(weighted.mean(x^2,invsigma))
  cat("m = ", Em, "+-", sigm, "\n")
  cat("c = ", Ec, "+-", sigc, "\n")
  cat("Coefficiente di correlazione = ", correl, "\n")

  if(!is.null(estX)){
    sb = sqrt( 1/(sum(invsigma)) )
    for (value in estX) {
      estY = Em * value + Ec
      sigy = sqrt( ((sb^2)/length(x)) * (1+ (value- meanx)^2/variance) )
      cat("\n", "Valore estrapolato:", "\n")
      cat("\t", "x: ", value, "\n")
      cat("\t", "y: ", estY, "+-", sigy, "\n")
    }
  }

  if(verbose){
    cat("\n", "> Media pesata delle x:", meanx , "\n")
    cat("> Media pesata delle y:", meany, "\n")
    cat("> Media pesata di xy:", weighted.mean(x*y,invsigma), "\n")
    cat("> Media pesata di x^2:", weighted.mean(x^2,invsigma), "\n")
    cat("> Varianza di x:", variance, "\n")
    cat("> Covarianza tra x e y:", weighted.mean(x*y, invsigma) - meanx*meany, "\n")
    cat("> Sommatoria 1/sigma^2:", sum(1/sigma^2),"\n")
  }


  plot(x, y)
  abline(Ec, Em)
  grid()


  invisible( list(x = x, y=y, m=Em, sigma.m=sigm, c=Ec, sigma.c=sigc, correl=correl) )
}