# data un valore di valori attesi, matrice di covarianza e
# vettori condizionanti ritorna vettori condizionati

norm_mult_cond <- function(mu, V, x.c, full=FALSE) {
  out <- NULL
  # Controlli: 
  n <- length(mu)
  # 1) dimensioni di V
  if ( sum(dim(V) != n) ) {
    cat( sprintf("dimensioni V incompatibile con lunghezza mu\n") )
    return(out)
  }
  # 2) V deve essere definita positiva (nessun autovettore <= 0)
  if( sum( eigen(V)$values <= 0 ) >  0) {
    cat( sprintf("V non è definita positiva\n") )
    return(out)
  }
  # numero di variabili condizionanti
  nc <- length(x.c[!is.na(x.c)])
  # casi particolari/anomali
  if( (length(x.c) > n) | (nc > n) ) {
    cat( sprintf("x.c ha più elementi di mu\n") )
    return(out)
  } else if (nc == 0) {  # non c'è alcuna condizione
    out$mu <- mu
    out$V  <- V
    return(out)
  } else if(nc == n) {
    out$mu <- x.c    # valori esatti
    out$V  <- NULL   # matrice di covarianza non ha senso
    return(out)
  }

  # calcoliamo vettore e matrice condizionate
  # Vedi  http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Bivariate_conditional_expectation
  v.c <- which(!is.na(x.c))  # variabili condizionanti
  v   <- which(is.na(x.c))   # variabili di interesse 
  V11 <- V[v, v]      # Sigma_11 di Wiki
  V22 <- V[v.c, v.c]  # Sigma_22  "  "
  V12 <- V[v, v.c]    # Sigma_12  "  "
  V21 <- V[v.c, v]    # Sigma_21  "  "
  mu.cond <- mu[v] + V12 %*% solve(V22) %*%  (x.c[!is.na(x.c)] - mu[v.c])
  V.cond  <- V11 - V12 %*% solve(V22) %*% V21
  if(!full) { # rimanda indietro solo la parte di interesse
    out$mu <- as.vector(mu.cond)
    out$V  <- V.cond
  } else {
    mu1 <- mu
    V1 <- V
    mu1[v]   <- mu.cond
    mu1[v.c] <- x.c[!is.na(x.c)]
    V1[v, v] <- V.cond
    V1[v.c, v.c] <- 0
    V1[v, v.c]   <- 0
    V1[v.c, v]   <- 0
    out$mu <- as.vector(mu1)
    out$V  <- V1
  }
  return(out)
}