function [distrib, valx, valy] = distrib_dos_va(xpart,xval,ypart,yval)
%  [distrib,valx,valy] = distrib_dos_va(xpart,xval,ypart,yval)
%  Dadas dos variables aleatorias en forma particional, hallar su distribucion
%  Las variables estan definidas sobre el espacio muestral [0,1],
%     con una distribucion uniforme
%  xpart: vector (fila o columna) de puntos que definen la particion de X
%         1er elemento: 0, ultimo elemento: 1, orden ascendiente
%  ypart: vector (fila o columna) de puntos que definen la particion de Y
%         1er elemento: 0, ultimo elemento: 1, orden ascendiente
%  xval:  valores para los distintos elementos de la particion de X
%  yval:  valores para los distintos elementos de la particion de Y
%  valx: valores _distintos_ que toma la variable aleatoria X (vect columna)
%  valy: valores _distintos_ que toma la variable aleatoria Y (vect fila)
%  distrib: distribucion conjunta del vector (X, Y), 
%    ie, probabilidad de cada par de valores de los vectores valx y valy
% 
%  Ejemplo:  xpart = [0 1/4 1/3 1/2 3/4 1]  xval =  [10  15  20  15  10]
%            ypart = [0 1/3 5/6 1]  yval = [3 2 3]
%  El output es: valx = [10; 15; 20], valy = [2 3]
%                distrib = [1/12 5/12; 3/12 1/12; 2/12 0]

%% Funcion de Matlab/Octave -- Ricard Torres, 2015 (sin acentos)

% En primer lugar, comprobamos que los argumentos sean correctos
if ( nargin ~= 4 )
   mostrar_uso() ;
   error('La funcion debe tener cuatro argumentos') ;
end

if ~ and( isvector(xpart), and(isvector(xval), and(isvector(ypart), ...
     isvector(yval) ) ) )
   mostrar_uso() ;
   error('Los argumentos de la funcion deben de ser vectores') ;
end

lpx = length(xpart); lvx = length(xval);
lpy = length(ypart); lvy = length(yval);

if ~ and( ( lvx == lpx - 1 ),  ( lvy == lpy - 1 ) ) 
   mostrar_uso() ;
   error('Debe de haber un valor para cada elemento de la particion') ;
end

if ~ and( (xpart(1) == 0), and( (xpart(end) == 1), and( (ypart(1) == 0), ...
    (ypart(end) == 1) ) ) ) 
   mostrar_uso() ;
   error('Cada particion debe empezar con 0 y acabar con 1') ;
end

is_ascending = true;
for i = 2 : lpx
    if xpart(i) <= xpart(i-1), is_ascending = false; end
end
for i = 2 : lpy
    if ypart(i) <= ypart(i-1), is_ascending = false; end
end

if ~ is_ascending
   mostrar_uso() ;
   error('El vector que define cada particion debe tener orden ascendiente') ;
end

% Ahora desarrollamos los calculos propiamente dichos
valx = unique(sort(xval));
valy = unique(sort(yval));
lx = length(valx);
ly = length(valy);
distrib = zeros(lx,ly);
% Asegurar que ambos vectores de particion son vectores fila
xpart = reshape(xpart,1,lpx);
ypart = reshape(ypart,1,lpy);
% Hallar la menor particion que refina tanto a xpart como a ypart
partcomun = unique(sort([xpart,ypart]));
% Hallar probabilidades conjuntas para los valores de ambas particiones
for i = 2 : length(partcomun)
    j = 2;
    while ( xpart(j) < partcomun(i) )
          j = j + 1;
    end
    k = 2;
    while ( ypart(k) < partcomun(i) )
          k = k + 1;
    end
    probconj(j-1,k-1) = partcomun(i) - partcomun(i-1);
end
% Hallar la distribucion de valores sin repeticiones
for i = 1 : lx
    for k = 1 : lvx
        if ( valx(i) ==  xval(k) )
           for j = 1 : ly
               for l = 1 : lvy
                   if ( valy(j) ==  yval(l) )
                      distrib(i,j) = distrib(i,j) + probconj(k,l);
                   end
               end
           end
        end
    end
end
% Dar valores de X como vector columna
valx = reshape(valx,lx,1);
% Dar valores de Y como vector fila
valy = reshape(valy,1,ly);

end

% Mensaje de ayuda en caso de error
function mostrar_uso()
fprintf('Uso: [distrib,valx,valy] = distrib_dos_va(xpart,xval,ypart,yval)\n');
fprintf('Dadas dos variables aleatorias en forma particional, ');
fprintf('hallar su distribucion\n');
fprintf('Las variables estan definidas sobre el espacio muestral [0,1],\n');
fprintf('  con una distribucion uniforme\n');
fprintf('xpart: vector (fila o columna) de puntos que definen la particion de X\n');
fprintf('    1er elemento: 0, ultimo elemento: 1, orden ascendiente\n');
fprintf('ypart: vector (fila o columna) de puntos que definen la particion de Y\n');
fprintf('    1er elemento: 0, ultimo elemento: 1, orden ascendiente\n');
fprintf('xval:  valores para los distintos elementos ');
fprintf(' de la particion de X\n');
fprintf('yval:  valores para los distintos elementos ');
fprintf(' de la particion de Y\n');
fprintf('valx: valores _distintos_ que toma la variable aleatoria X ');
fprintf('(vect columna)\n');
fprintf('valy: valores _distintos_ que toma la variable aleatoria Y ');
fprintf('(vect fila)\n');
fprintf('distrib: distribucion conjunta del vector (X, Y)\n');
fprintf(' ie, probabilidad de cada par de valores de los vectores valx y valy\n');
fprintf('\nEjemplo:  xpart = [0 1/4 1/3 1/2 3/4 1] xval =  [10  15  20  15  10]\n');
fprintf('            ypart = [0 1/3 5/6 1]  yval = [3 2 3]\n');
fprintf('El output es: valx = [10; 15; 20], valy = [2 3]\n');
fprintf('           distrib = [1/12 5/12; 3/12 1/12; 2/12 0]\n');
end