function grafica_dos_lorenz(na, xa, nb, xb)
%% function grafica_dos_lorenz(na, xa, nb, xb)
%% Graficar las curvas de Lorenz de dos distribuciones de ingresos finitas
%% ni: son las frecuencias absolutas (números) de cada categoria
%% xi: son las categorias de ingreso (en estricto orden ascendiente)
%%  Ricard Torres <rtorres@itam.mx>

   la = length(xa);
   lb = length(xb);
   if la ~= length(na) || lb ~= length(nb)
      error('Categorias y frecuencias deben tener la misma longitud');
   end
   if or( any( na <= 0 ), or( any( xa <= 0 ), ...
             or( any ( nb <= 0 ), any ( xb <= 0 ) )))
      error('Categorias y frecuencias deben ser estrictamente positivos');
   end
   if any(xa ~= sort(unique(xa))) || any(xb ~= sort(unique(xb)))
      error('Categorias deben estar en orden ascendiente estricto');
   end
   [Fa, Sa] = lorenz_disc(na,xa); % Curva de Lorenz
   [Fb, Sb] = lorenz_disc(nb,xb);
   diag = [0, 1]; % diagonal
   figure; % crear nueva figura
   plot(diag,diag,'--','Color','Green','LineWidth',1);
   hold on;
   ha = plot(Fa,Sa,'-','Color','Blue','LineWidth',2);
   plot(Fa,Sa,'o','Color','Blue','LineWidth',2); % Marcar puntos
   hb = plot(Fb,Sb,'-','Color','Red','LineWidth',2);
   plot(Fb,Sb,'o','Color','Red','LineWidth',2); % Marcar puntos
   axis([-0.05,1.05,-0.05,1.05]);
   xlabel('Frecuencias acumuladas poblacion'); 
   ylabel('Proporciones acumuladas ingreso');
   title('Curvas de Lorenz');
   legend([ha,hb],{'Distrib A','Distrib B'},'Location','NorthWest');
   hold off
end