# Script de Julia (http://julialang.org/) para calcular 
# varianzas y covarianza de un vector aleatorio (X,Y) y de su suma
# Ricard Torres <rtorres@itam.mx>
# Para generar datos, cargamos la libreria de distribuciones
# (eso toma su tiempo, pues de otro modo la ejecucion seria instantanea)
using Distributions
du = DiscreteUniform # Usar dist uniforme discreta para generar datos
# Generar valores v.a. y ponerlos en orden ascendiente 
X = sort(unique(rand(du(100,200),6))) 
Y = sort(unique(rand(du(50,150),8))) 
lx = length(X) # Longitud vector X (puede ser menor)
ly = length(Y)
# Matriz de la distribucion conjunta del vector (X,Y)
P = rand(du(0,100),lx,ly)  # frecuencias absolutas
P = (1/sum(P)) * P        # frecuencias relativas (probabilidades)
# Distribuciones marginales (hay que convertirlas a vectores)
PX = vec(sum(P,2)) 
PY = vec(sum(P,1)) 
# Valores esperados (medias)
EX = dot(PX, X) 
EY = dot(PY, Y) 
# Variables aleatorias centradas alrededor de la media
XC = X - EX * ones(size(X))
YC = Y - EY * ones(size(Y))
# Varianzas
VX = dot( PX, XC.^2 )
VY = dot( PY, YC.^2 )
# Covarianza del vector (X,Y)
# (aseguramos compatibilidad de dimensiones mediante 'reshape')
CovXY = reshape(XC,1,length(XC)) * P * reshape(YC,length(YC),1)
# Truco para convertir CovXY de vector de longitud 1 a escalar (julia quirk)
CovXY = sum(CovXY) 
# Coeficiente de correlacion
CorrXY = CovXY / sqrt(VX * VY)
# Ahora calcularemos la suma, su distribucion y su varianza directamente
# En primer lugar, creamos una matriz con las dimensiones de P que
#  contiene las respectivas sumas
SS = reshape(X,lx,1) * ones(1,ly) + ones(lx,1) * reshape(Y,1,ly)
# Ahora pasamos las sumas y su distribucion a vectores
S = reshape(SS,prod(size(SS)))
PS = reshape(P,prod(size(P)))
# Esperanza de la suma
ES = dot( PS, S)
# Suma centrada alrededor de la media
SC = S - ES * ones(size(S))
# Varianza de la suma
VS = dot( PS, SC.^2 )
# Se debe cumplir que VS = VX + VY + 2 * CovXY, vamos a comprobarlo
tol = 1e-10
if abs( VX + VY + 2 * CovXY - VS ) > tol
   println("Ha habido algun error en los calculos de VS")
else
   println("El calculo de la varianza de la suma coincide con los dos metodos")
   println("X = ", X)
   println("Y = ", Y)
   println("P = ") 
   print(round(P,5))  # mostrar solamente 5 decimales
   println("")
   println("Var(X) = ", VX)
   println("Var(Y) = ", VY)
   println("Cov(X,Y) = ", CovXY)
   println("Corr(X,Y) = ", CorrXY)
   println("Var(X+Y) = ", VS)
   println("Var(X) + Var(Y) + 2 * Cov(X,Y) = ", VX+VY+2*CovXY)
end