# propagazione di matrice di covarianze mediante derivate
 
# genera la matrice delle espressioni delle derivate
mat.der <- function (Y, X) {
    mat <- expression()
    for (i in 1:length(X)) {                      # l'ordine e' per avere
       for (j in 1:length(Y)) {                   # in uscita righe in Y
          mat <- c(mat,  D(Y[[j]],names(X)[i]))
       }
    }
    dim(mat) <- c(length(Y),length(X))            # vedi sopra...
    return(mat)
}   
 
# facciamo la funzione che calcola i valori numerici
mat.der.num <- function (mat, X) {
    attach(X)                    # rende visibile gli elementi di X,
                                 # ovvero a=... e b=...
    matn <- numeric()            # matrice vuota 
    for  (i in 1:length(mat)) {
      matn[i] <- eval(mat[i])    # si valutano tutte derivate
    }
    dim(matn) <- dim(mat)        # assegnando le dimensioni opportune
    return(matn)
}   
 
# 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) ) )
}
 
# valori attesi (ovviamente in approx lineare)
E.Y <- function(Y, X) {
    attach(X)
    ey <- numeric()
    for (i in 1:length(Y))  ey[i] <- eval(Y[[i]]) 
    return(ey)
}
 
# 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)
 
# Variabili di output
ExpY <- E.Y(Y, X)
mat  <- mat.der(Y, X)
matn <- mat.der.num(mat, X)
covY <- matn %*% ( covX  %*% t(matn) )
 
sigmaY <- sqrt( diag(covY) )
corY   <- covY / sqrt( diag(covY) %*% t(diag(covY)) )