%% Variables aleatorias multidimensionales (ejemplo Matlab)
% © 2020, Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo en: http://personales.upv.es/asala/YT/V/vamdm.html
% Este código ejecutó correctamente en Matlab R2019b
% Table of Contents
% Generación de datos
% Estudio de las distribuciones marginales
%     Media y Varianza
%     Momentos superiores marginales
% Covarianza y correlación entre las variables
%
%% Generación de datos

N=1000;
[x1,x2,x3]=formulasecreta(N)

%% Estudio de las distribuciones marginales

%representemos graficamente los histogramas "marginales"
subplot(1,3,1)
histogram(x1), grid on
subplot(1,3,2)
histogram(x2), grid on
subplot(1,3,3)
histogram(x3), grid on

%% Media y Varianza

%calculemos las medias marginales:
media1=sum(x1)/N
media2=sum(x2)/N
media3=sum(x3)/N

%calculemos las variables incrementales:
x1inc=x1-media1;
x2inc=x2-media2;
x3inc=x3-media3;

%las varianzas (muestrales) marginales:
var1=sum(x1inc.^2)/(N-1)
var2=sum(x2inc.^2)/(N-1)
var3=sum(x3inc.^2)/(N-1)
var(x1) %también lo hace Matlab, claro!.
%desviación típica
s1=std(x1)
s2=std(x2)
s3=std(x3)
%en distribución normal, el 98% de las muestras están entre +/- 2.33 veces
%su desv. típica
lll=-norminv(0.02/2)
length(find(-lll*s1<=x1inc & x1inc<=lll*s1))/N %tiene más outliers "extremos"...
length(find(-lll*s2<=x2inc & x2inc<=lll*s2))/N
length(find(-lll*s3<=x3inc & x3inc<=lll*s3))/N

%% Momentos superiores marginales

% skewness (sesgo, asimetría, mediaxstd^3), kurtosis (prob. valores extremos, media xstd^4)
% distrib. normal tiene sesgo=0, curtosis=3;
skew=skewness([x1' x2' x3'])
kur=kurtosis([x1' x2' x3'])
%la variable 1 tiene mucha más curtosis que la distr. normal.

%% Covarianza y correlación entre las variables

%la covarianza entre x1 y x2:
cov12=sum(x1inc.*x2inc)/(N-1)
% mediante la matriz de datos:
X=[x1inc;x2inc;x3inc];
Sigma=X*X'/(N-1)
cov(X') %Matlab lo tiene implementado... muestras en "filas"
eig(Sigma)
%si normalizamos todas las vbles a media 0, varianza 1:
x1incstd=x1inc/sqrt(var1);
x2incstd=x2inc/sqrt(var2);
x3incstd=x3inc/sqrt(var3);
Xstd=[x1incstd;x2incstd;x3incstd];
X2=zscore(X')'; %Matlab también lo tiene preparado.
%estodacero=X2-Xstd %esto da CERO
%obtenemos la "pairwise-correlation" matrix con variables normalizadas
Sigmastd=cov(Xstd')
%dibujemos el pairwise scatterplot:
for i=1:3
    subplot(3,3,1+(i-1)*4)
    histogram(Xstd(i,:)), xlabel('x'+string(i)), grid on
    for j=(i+1):3
        subplot(3,3,(i-1)*3+j)
        plot(Xstd(j,:),Xstd(i,:),'.'), xlabel('x'+string(j)),  ylabel('x'+string(i))
        hold on
        plot(Xstd(j,:),Xstd(j,:)*Sigmastd(i,j),'-.','LineWidth',2)
        grid on
        hold off
    end
end

%%









%

function [x1,x2,x3]=formulasecreta(Nmuestras)
    p=randn(3,Nmuestras);
    q=[1 2 3;-2 -3 1;3 6 7]*p;
    x1=q(1,:).^3+randn(1,Nmuestras)*0.3;
    x2=q(2,:)*1.5;
    x3=q(3,:);
end