# misura lati interni di un triangolo con vincolo
#
#  G. D'Agostini, 23/03/2010
#-------------------------------------------------


mu.X    <- c(58, 73, 54)
Sigma.X <- matrix(rep(0,9), 3, 3)
Sigma.X[1,1] <-  2.0^2
Sigma.X[2,2] <-  2.0^2
Sigma.X[3,3] <-  2.0^2

C.m <- rbind( c(1,0,0), c(0,1,0), c(0,0,1), c(1,1,1) )

mu.Y <- C.m %*% mu.X
Sigma.Y <- C.m %*% Sigma.X %*% t(C.m)
sigma.Y <- sqrt(c(Sigma.Y[1,1],Sigma.Y[2,2],Sigma.Y[3,3],Sigma.Y[4,4]))
rho.Y <- Sigma.Y / ( sigma.Y  %*% t(sigma.Y))
# usa package mnormt (multidimensionale!)
#library(mnormt)

somma = 180.
dif <- as.vector(somma-mu.Y[4])
# parametri della distribuzione condizionata:
# http://en.wikipedia.org/wiki/Multivariate_normal_distribution

Sigma.Y.11 = Sigma.X

Sigma.Y.12 = matrix(rep(0,3), 3, 1)
Sigma.Y.12[1:3,1] = Sigma.Y[1:3,4]

Sigma.Y.21 = matrix(rep(0,3), 1, 3)
Sigma.Y.21[1,1:3] = Sigma.Y[4,1:3]

Sigma.Y.22 = matrix(Sigma.Y[4,4] , 1, 1)

mu.X.cond <- as.vector(mu.X +  Sigma.Y.12 %*% solve(Sigma.Y.22) %*% dif)
Sigma.X.cond <- Sigma.Y.11 -  Sigma.Y.12 %*% solve(Sigma.Y.22) %*% Sigma.Y.21

sigma.X.cond <- c(sqrt(Sigma.X.cond[1,1]),
                  sqrt(Sigma.X.cond[2,2]),
                  sqrt(Sigma.X.cond[3,3]) )

corr.X.cond <- Sigma.X.cond / (sigma.X.cond %*% t(sigma.X.cond))

cat(sprintf("\n Valori attesi degli angoli ricondizionati dalla somma: \n"))
print(mu.X.cond)
cat(sprintf("\n Deviazioni standard:\n"))
print(sigma.X.cond)
cat(sprintf("\n Matrice di covarianza:\n"))
print(Sigma.X.cond)
cat(sprintf("\n Matrice di correlazione:\n"))
print(corr.X.cond)
# da confrontare con i risultati di JAGS/rjags