  function [goodsamples,badsamples]=RejectionSampling(N,pfun,qfun,M,draw,range)
  
  % RejectionSampling      Create samples with the Rejection sampling
  %                        method
  %
  % [goodsamples,rejecedsamples]=RejectionSampling(N,pfun,qfun,draw,range)
  % N=# of accepted samples to create
  % pfun: Target original function (see define function)
  % qfun: Proposal function (see define function)
  % M: M factor
  % draw: 1, draw result information, 0 just sample
  % range: range values to draw
  % goodsamples: accepted samples
  % badsamples: rejected samples
  
  % APRENDIZAJE AUTOMÁTICO (G. Pajares y J. M. de la Cruz, Eds.)
  % Editorial: RA-MA
  % Capítulo 7 y 8: Monte Carlo I y II
  % Autores: Eva Besada-Portas, Sergey M. Plis

  
  
% ----------------------------------------------------------- SAMPLING
  
  remaining=N;
  goodsamples=[];goodchk=[];
  badsamples=[];badchk=[];
  
  while (remaining) %There are still elements to sample

    samples=qfun.fun(1,remaining,qfun.param{:}); %Sample from q
    q= qfun.fun(0,samples,qfun.param{:});        %q(samples)   
    p=pfun.fun(0,samples,pfun.param{:});             %p(samples)
  
    chk=M*rand(1,remaining).*q;                 %M*u*q
    good=  find(chk-p<0);                       %M*u*q<p
    bad=setdiff([1:remaining]',good);           %The others
    badsamples = [badsamples,samples(bad)];     %rejected samples
    badchk=[badchk,chk(bad)];                   %M*u*q of the rejected samples
    goodsamples=[goodsamples,samples(good)];    %accepted samples
    goodchk=[goodchk,chk(good)];                %M*u*q of the accepted samples
    remaining=remaining-length(good);           %How many samples need to be created to reach N
    %badchk     = badchk + pfun(badsamples, means,vars);
  end
  
  
% ----------------------------------------------------------- PLOTTING
  if ~draw
      return;
  end
  
  %Show the distribution, and a part of the good and bad samples
  subplot(311);
  
  group=min(100,N); %We only represent a subset because N can be too big
  
  %Represent part of the rejected samples
  hmb=length(badsamples);
  hmbd=floor(hmb/N*group); 
  ind=randperm(hmb);
  ind=ind(1:hmbd);
  plot(badsamples(ind), badchk(ind),'bs','MarkerFaceColor',[0 0 0.5],'MarkerSize',4);
  
  %Represent part of the accepted samples
  hold on
  ind=randperm(N);
  ind=ind(1:group);
  plot(goodsamples(ind), goodchk(ind),'*','Color',[0 0.5 0]);
  
  %Show the distributions
  r = pfun.fun(0,range,pfun.param{:}); 
  plot(range,r,'g','LineWidth',3);
  r = M*qfun.fun(0,range,qfun.param{:}); 
  plot(range,r,'r', 'LineWidth',2);
  h=legend('Rech.','Acept.','p(C)','M*q(C)');
  set(h,'FontSize',12);
  title(sprintf('a) M=%2.1f, Aceptadas=%d, Rechazadas=%d',M,N,length(badsamples)),'FontSize',12);
  %set(gca,'xtick',[])
  ylims = ylim();
  ylim([-0.05 ylims(2)]);
  xlim([min(range),max(range)]);
  hold off

  %Show the number of rejected and accepted samples
  subplot(312);
  HistogramCount(goodsamples,badsamples);
  h=legend('Acept.','Rech.');
  set(h,'FontSize',12);
  title('b) Numero de muestras','FontSize',12);
  xlim([min(range),max(range)]);
  
  %Show the weighted histogram
  subplot(313);
  HistogramWeighted(goodsamples,ones(1,N)/N);
  hold on;
  r= pfun.fun(0,range,pfun.param{:}); 
  plot(range,r,'g','LineWidth',3);
  h=legend('Hist.','p(C)');
  set(h,'FontSize',12);
  title('c) Histograma','FontSize',12);

  ylims = ylim();
  ylim([-0.05 ylims(2)]);  
  xlim([min(range),max(range)]); 
  xlabel('Concentracion (%)','FontSize',12);

  