#----------------------------------------------------------
# covariance matrix and correlation matrix of sum and difference
# (exercise on the measurements of the sides of the A4 foil)
#
#  GdA 26/2/2016
#------------------------------

m.a = 29.73
m.b = 21.41

sd.a = 0.03    # then change to 0.04 and to 0.05 to see the difference!
sd.b = 0.04

rho.ab = 0.   # then change to +0.8 to see the difference

cat("\nInput  ---------------------------------------------------------------- \n")
cat(sprintf("a = (%.2f +- %.2f) cm\n", m.a, sd.a ))
cat(sprintf("b = (%.2f +- %.2f) cm\n", m.b, sd.b ))
cat(sprintf("rho(a,b) = %.2f \n", rho.ab ))

#-------------------------------------------------------------------------------
# 1) just using the formulae  (coefficients highlighted) 

# standard deviations,covariance and correlation coefficient
sd.s = sqrt( sd.a^2 + sd.b^2 + 2*(1)*(1)* (rho.ab*sd.a*sd.b) ) 
sd.d = sqrt(sd.a^2 + sd.b^2 + 2*(1)*(-1)* (rho.ab*sd.a*sd.b) ) 
cat("\nPropagation ---------------------------------------------------------- \n")
cat(sprintf("s = (%.2f +- %.2f) cm\n", m.a+m.b, sd.s ))
cat(sprintf("d = (%.2f +- %.2f) cm\n", m.a-m.b, sd.d ))

# The coefficients have been made explicit!
cov.sd = 1*1*sd.a^2 + 1*(-1)*sd.b^2 +
         ( 1*(-1) + 1*1 ) * ( rho.ab * sd.a * sd.b )   # IDENTICALLY NULL !!   
                                                      # correlations are tricky!
rho.sd = cov.sd / (  sd.s * sd.d )         # rho.ab DOES NOT influence cov.sd,
                                           # but it scales rho.sd 
                                           # because it DOES affect variances!
cat(sprintf("cov(s,d) = %.3e cm^2\n", cov.sd ))
cat(sprintf("rho(s,d) = %.3f \n", rho.sd))


#---------------------------------------------------------------------------------
# 2) Brute force, by MC   (see bivariate_gaussian_r.R for bivariate Gaussian)
n = 20000
a <- rnorm(n, m.a, sd.a)
b <- rnorm(n, m.b + rho.ab*sd.b/sd.a*(a-m.a), sd.b*sqrt(1-rho.ab^2))
s <- a+b
d <- a-b

cat("\nBrute force (direct sampling MC) -------------------------------------\n")
cat(sprintf("s = (%.2f +- %.2f) cm\n", mean(s), sd(s) ))
cat(sprintf("d = (%.2f +- %.2f) cm\n", mean(d), sd(d) ))
cat(sprintf("cov(s,d) = %.3e cm^2\n", cov(s,d) ))
cat(sprintf("rho(s,d) = %.3f \n", cor(s,d)))
plot(s,d, pch=19, cex=0.1, col='blue',
     main=sprintf("rho(a,b) = %.2f",cor(s,d)))

#------------------------------------------------------------------------------
# 3) Use compact matrice method

# transformation matrix
C <- rbind(c(1,1),c(1,-1))          

# Covariance matrix of the input quantities (a and b)
V <- matrix(rep(0,4),c(2,2))            # fill a 2x2 matrix with zeros
diag(V) = c(sd.a^2, sd.b^2)             # fill the diagonal termes with variances
V[1,2] = V[2,1] = rho.ab * sd.a * sd.b  # fill the offdiagonal with covariances

cat("\nPropagation by covariance matrix -------------------------------------\n")
cat("transformation matrix C\n")
print(C)
cat("\n'Input' covariance matrix V\n")
print(V)
cat("\n'Input' correlation matrix V\n")
print( V / outer( c(sd.a,sd.b), c(sd.a,sd.b) ) )


Vy <- C %*% V %*% t(C)
cat("\n\n'Output' covariance matrix V\n")
print(Vy)
cat("\n'Output' standard deviations V\n")
sd.y <- sqrt(diag(Vy))               # note how "diag()" is bidirectional!
print(sd.y)
cat("\n'Output' correlation matrix matrix V\n")
Cor.y <- Vy / outer(sd.y, sd.y)
print(Cor.y)