% Runs all the Auxiliar Samplers once

% 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


close all;%clear all;

%--------------------------------PARAMETROS------------------------------
 N       = 10000;          %Number of samples
 proposal=0;               % Selected qfun
                           % 0: uniform
                           % 1: gaussian
 draw=1;                   % Draw results inside each method
 range=0:0.1:100;          % Range values for values. 

 pfun=DefineFunction('pC'); %pfun funtion 
  
 switch(proposal)
      case 0,               %Uniform qfun in the limit range
        limit=[10,90];
        qfun=DefineFunction('quniform',limit);
        M=7.5;              %M for rejection sampling
        
      case 1,               %Gaussian qfun with mu mean and sigma std
        mu    = 43.5;
        sigma = 15;
        qfun=DefineFunction('qgaussian',mu,sigma);
        M = 3.53;           %M for rejection sampling
 end
  
  
  %Draw the original distribution
  figure('Name','Inicial');
  r = pfun.fun(0,range,pfun.param{:});
  plot(range,r,'g','LineWidth',3);
  xlabel('Concentracion (%)','FontSize',12);
  ylabel('p(Concentracion)','FontSize',12);
  title('Funcion de probabilidad original','FontSize',12);

  
  %Rejection sampling
  figure('Name','Rejected Sampling');
  [samplesRS,badsamples]=RejectionSampling(N,pfun,qfun,M,draw,range);
 
  %Importance Sampling
  figure('Name','Importance Sampling');
 [samplesIS,weightsIS]=ImportanceSampling(N,pfun,qfun,draw,range);
 
 %Markov Chaing 
 figure('Name','Markov Chain');
 s0=rand*10+80;
 [samplesMC]=MarkovChain(N,DefineFunction('ChainIndependentSampler',pfun,qfun),s0,draw,range,pfun);
 
 figure('Name','Markov Chain');
 [samplesMC]=MarkovChain(N,DefineFunction('ChainMetropolis',pfun,DefineFunction('qgaussian',5)),s0,draw,range,pfun);
 

 %Resampling with weights
 figure('Name','Reweight with weights');
 resample=weightsIS;
 [samplesRE,weightsRE]=WeightedResampling(samplesIS,weightsIS,resample,draw,range,pfun);
 
 
 figure('Name','Reweigth with uniform');
 resample=quniform(0,samplesIS,[10 90]);
 [samplesRE,weightsRE]=WeightedResampling(samplesIS,weightsIS,resample,draw,range,pfun);
 
 
 figure('Name','Reweight with gaussian');
 resample=qgaussian(0,samplesIS,43.5,15);
 [samplesRE,weightsRE]=WeightedResampling(samplesIS,weightsIS,resample,draw,range,pfun);
 
 