#----------------------------------------------
# linear fit done with JAGS repeated with Least Squares
# using matrix formalism
#
# GdA 27 marzo 2015
#----------------------------------------------
# function pausa
pausa <- function() {
  cat ("\n >> Guarda il plot e dai enter per continuare\n")
  scan()
}

data <- source("fit_data_OB.txt")$value
attach(data, warn=FALSE)  # make the elements of the list data 'visible'

plot(x, y, xlim=c(0, max(x)), ylim=c(-5, max(y)), col='blue')

# use lm()
fit.lm <-lm(y ~ x)

abline(fit.lm, lty=2, col='blue')

# -- use least square formulae -------------------------

# design matrix
X <- matrix(rep(1,n),c(n,1))
X <- cbind(X, x)

cm <- as.vector( solve(t(X)%*%X) %*% t(X) %*% y )   # best estimates ("beta")
abline(cm, lty=3, col='red')

# residuals
r <- y - X %*% cm
sigma <- as.vector(sqrt( t(r)%*%r / n))

# covariance matrix
Cov <- sigma^2 * solve( t(X) %*% X)
sigma.c <- sqrt(Cov[1,1])
sigma.m <- sqrt(Cov[2,2])
cor.cm <- Cov[1,2] / (sigma.c*sigma.m)