# propagazione di matrice di covarianze mediante MC
# vedi prop_cov.R per versione mediate formule standard di propagazione
# (mediante le derivate)
 
# libreria che contiene la normale multivariata
library(mvtnorm)
 
# Matrice di covarianza da sigma e correlazioni.
# In principio inutile (una sola riga di codice), ma fa anche dei controlli
build.cov <- function(corX, sigmaX) {
    if ( nrow(corX) != ncol(corX) )            return() # quadrata?
    if ( length(sigmaX) != nrow(corX) )        return() # dimens. compatibili?
    if ( any ( (t(corX) == corX) == FALSE ) )  return() # simmetrica?
    if ( any( corX > 1) ||  any( corX < -1) )  return() # valori ragionevoli?
    if ( any( diag(corX) != 1 ) )              return() # diagonale unitaria? 
    return ( corX * (sigmaX %*% t(sigmaX) ) )
}
 
# esempio di controllo
# espressioni delle variabili di 'output'
Y <- list( expression( 2*(a+b) ),         # perimetro
          expression( a*b ),              # area
          expression( sqrt(a^2+b^2) ) ,   # diagonale
          expression(a/b)                 # rapporto lati (ovvero tangente)
          )
 
# variabili di input
# valori attesi, sigma, matrice di correlezione e matrice di covarianza
X <- list(a=29.73, b=21.45)
sigmaX <- c(0.03, 0.04)
# inizialmente di correlazione matrice diagonale
corX <- matrix( rep(0, length(X)^2), nrow=length(X) )
diag(corX) <- 1
# eventualmente possiamo aggiungere delle correlazioni, ad es,
corX[1,2] <- corX[2,1] <- 0.0    # sostituire con valori [-1, +1]
# matrice di covarianza
covX <- build.cov(corX, sigmaX)
 
#--------
n <- 10000
x <- rmvnorm(n, c(X$a, X$b), covX)
a <- x[,1]
b <- x[,2]
 
y <- c( eval(Y[[1]]), eval(Y[[2]]), eval(Y[[3]]), eval(Y[[4]]) )
dim(y) <- c(n, length(Y))
 
ym <- numeric()
yv <- matrix(rep(0, length(Y)^2), nrow=length(Y))
yr <- matrix(rep(0, length(Y)^2), nrow=length(Y))
for (i in 1:length(Y) ) {
  ym[i] <- mean(y[,i])                # medie
  for (j in 1:length(Y)) {
    yv[i,j] <- cov(y[,i],y[,j])       # matrice di covarianza
  }
  ys <- sqrt(diag(yv))                # sigma
  yr <- yv / (ys %*% t(ys) )          # matrice di correlazione
}