function muestra = inversa_distrib_discreta(num_aleat,probs,vals)
%  muestra = inversa_distrib_discreta(num_aleat,probs,vals)
%  Dada una distribucion discreta, y un numero (vector) aleatorio 
%  entre 0 y 1, devolver el cuantil (vector) correspondiente
%  Si 'num_aleat' es una muestra de una distribucion uniforme en [0,1],
%  'muestra' es una muestra del mismo tama~no de la distribucion discreta
%  num_aleat: muestra aleatoria de la distrib uniforme en [0,1]
%  vals:  valores que toma la distribucion discreta en cuestion
%  probs: probabilidades correspondientes a los valores anteriores
%  muestra: muestra aleatoria de la distribucion discreta
% 
%  Ejemplo:  probs = [ 1/4 1/3 1/12 1/6 1/6 ]
%            vals   =  [  10  15  20   30  40 ]
%            num_aleat = [ 0.25 0.87 ]
%  El output es: muestra = [ 10 40 ] 

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

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

if ~ and( isvector(num_aleat), and( isvector(probs), isvector(vals) ))
   mostrar_uso();
   error('Los argumentos de la funcion deben de ser vectores');
end

ln = length(num_aleat); 
lp = length(probs); 
lv = length(vals);

if ( lv ~= lp )
   mostrar_uso();
   error('Debe haber el mismo numero de probabilidades y valores') ;
end

if any( probs < 0 )
   mostrar_uso() ;
   error('Las probabilidades deben ser no negativas') ;
end

tolerancia = 1e-10;
if abs( sum(probs) - 1 ) > tolerancia
   mostrar_uso() ;
   error('Las probabilidades deben sumar 1') ;
end

if ~ all( num_aleat >= 0 )
   mostrar_uso() ;
   error('Los numeros aleatorios deben ser no negativos') ;
end

if ~ all( num_aleat <= 1 )
   mostrar_uso() ;
   error('Los numeros aleatorios no pueden exceder 1') ;
end

% Ahora desarrollamos los calculos propiamente dichos

% Escribir 'vals' en orden ascendiente si no lo esta, 
% cambiando 'probs' para que se correspondan
if any( vals ~= sort(vals) )
   for i = 1 : lv
       [provis , j] = min(vals(i:lv)); % solo usamos j, no provis
       if ( i ~= j )
          provis = vals(i); % intercambiar vals(i) y vals(j)
          vals(i) = vals(j);
          vals(j) = provis;
          provis = probs(i); % intercambiar probs(i) y probs(j)
          probs(i) = probs(j);
          probs(j) = provis;
       end
   end
end
% Eliminar repeticiones en 'vals' si las hay, 
% cambiando 'probs' para que se correspondan
if ( length(unique(vals)) < lv ) 
   i = 1;
   j = 1;
   while ( i <= lv )
       nv(j) = vals(i);
       np(j) = probs(i);
       while ( ( i < lv ) && ( vals(i) == vals(i + 1) ) ) 
          np(j) = np(j) + probs(i+1) ;
          i = i + 1;
       end
       i = i + 1;
       j = j + 1;
   end
   vals = nv;
   probs = np;
   lv = length(vals);
end

% Funcion de distribucion
fdist = cumsum(probs);

% Aplicamos la inversa generalizada de la funcion de distribucion:
% si fdist(j-1) < num_aleat <= fdist(j), entonces muestra = vals(j)
for i = 1 : ln   % para cada numero aleatorio, hallar su imagen inversa
    for j = 1 : lv
        if ( num_aleat(i) <= fdist(j) )
           muestra(i) = vals(j);
           break;
        end
     end
end
        
end

% Mensaje de ayuda en caso de error
function mostrar_uso()
fprintf('Uso: muestra = inversa_distrib_discreta(num_aleat,probs,vals)\n');
fprintf('Dada una distribucion discreta, y un numero (vector) aleatorio\n');
fprintf('entre 0 y 1, devolver el cuantil (vector) correspondiente\n');
fprintf('Si ''num_aleat'' es una muestra de una distribucion uniforme en [0,1],\n');
fprintf('''muestra'' es una muestra del mismo tama~no de la distribucion discreta\n');
fprintf('num_aleat: muestra aleatoria de la distrib uniforme en [0,1]\n');
fprintf('vals:  valores que toma la distribucion discreta en cuestion\n');
fprintf('probs: probabilidades correspondientes a los valores anteriores\n');
fprintf('muestra: muestra aleatoria de la distribucion discreta\n');
fprintf('\nEjemplo:  probs = [ 1/4 1/3 1/12 1/6 1/6 ]\n');
fprintf('             vals =  [  10  15  20   30  40 ]\n');
fprintf('             num_aleat = [ 0.25 0.87 ]');
fprintf('El output es: muestra = [ 10 40 ]\n');
end