function samples=MarkovChain(N,chainfun,samples,draw,range,pfun)

  % Markov chainfun method
  %
  % samples=MarkovChain(N,chainfun,s0,draw,range)
  % N=# of samples
  % chainfun: chain function (see define function)
  % s0: initial sample
  % draw: 1, draw result information, 0, just sample
  % range: range values to draw
  % pfun: q(C) to draw
  % samples: accepted 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
 
  i=2;
  for j=2:N,
      samples(i)=chainfun.fun(samples(i-1),chainfun.param{:});
      if j>1 i=i+1; 
      else
          samples(i-1)=samples(i);
      end
  end
  
 
% ----------------------------------------------------------- PLOTTING
 if ~draw
      return;
 end
  
 %Draw part of the chainfun (last groups of experiments)
  subplot(311)
  hmdraw=min(N,5000);
  drawing=N-hmdraw+1; %We only plot last 500 examples
  plot(samples(drawing:N),N:-1:drawing,'o-','Color',[0.5 .5 .5],'MarkerEdgeColor',[0 0.4 0.4],'MarkerFaceColor',[0 0.8 0.8]);
  title(sprintf('a) Ultimas %d muestras',hmdraw),'FontSize',12);
  set(gca,'ytick',[])
  xlim([min(range),max(range)]);
  
  %Histogram count
  subplot(312);
  HistogramCount(samples);
  title('b) Numero de muestras','FontSize',12);
  xlim([min(range),max(range)]);
  
  %Weighted histogram and original distribution
  subplot(313);
  HistogramWeighted(samples,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);
  hold off
  