% 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


clear all;
Data;

%Representing the p(C,T)
CRange=1:100;
TRange=1:1500;
for C=CRange
    for T=TRange
        if T<s1
            P(T,C)=1/sqrt(2*pi*c1)*exp(-(C-m1)^2/2/c1)/s5;
        elseif T<s2
            P(T,C)=1/sqrt(2*pi*c2)*exp(-(C-(m1+(m2-m1)/(s2-s1)*(T-s1)))^2/2/c2)/s5;           
        elseif T<s3
            P(T,C)=1/sqrt(2*pi*c3)*exp(-(C-m2)^2/2/c3)/s5; 
        elseif T<s4
            P(T,C)=1/sqrt(2*pi*c4)*exp(-(C-(m2+(m3-m2)/(s4-s3)*(T-s3)))^2/2/c4)/s5;           
        elseif T<=s5        
            P(T,C)=1/sqrt(2*pi*c5)*exp(-(C-m3)^2/2/c5)/s5; 
        end
    end
end
figure;
[CC,TT]=meshgrid(CRange,TRange);
mesh(CC,TT,P);
xlabel('Concentracion (%)','FontSize',12);
ylabel('Tiempo (horas)','FontSize',12);
title('p(C,T)','FontSize',12)
view(-57.5,66);

        
%Generating the example for a single experiment, following the
%distributions
h1=gsamp(m1,c1,s1)'; %First part
h2=gsamp(0,c2,s2-s1)'+linspace(m1,m2,s2-s1); %Second part
h3=gsamp(m2,c3,s3-s2)'; %Third part
h4=gsamp(0,c4,s4-s3)'+linspace(m2,m3,s4-s3); %Fourth part
h5=gsamp(m3,c5,s5-s4)'; %Fifth part
t=0:s5-1;
h=[h1,h2,h3,h4,h5];  %values at the different time steps
figure
plot(t,h,'LineWidth',2);
h=title('Concetracion en un experimento','FontSize',12);
xlabel('Tiempo (horas)','FontSize',12);
ylabel('Concentración (%)','FontSize',12);
ylim([0,100])

%Mean Concentration at each time step
h1=linspace(m1,m1,s1);
h2=linspace(m1,m2,s2-s1);
h3=linspace(m2,m2,s3-s2);
h4=linspace(m2,m3,s4-s3);
h5=linspace(m3,m3,s5-s4);
h=[h1,h2,h3,h4,h5];
figure
plot(t,h,'LineWidth',2,'Color','r');
h=title('Concetracion media','FontSize',12);
xlabel('Tiempo (horas)','FontSize',12);
ylabel('Concentración media (%)','FontSize',12);
ylim([0,100])
            
%P(C) Marginal
%P(C)=int(P(C,T))
X=0:.1:100;
P1=s1/sqrt(2*pi*c1)*(exp(-(X-m1).^2/2/c1))/s5;
P3=(s3-s2)/sqrt(2*pi*c3)*(exp(-(X-m2).^2/2/c3))/s5;
P5=(s5-s4)/sqrt(2*pi*c5)*(exp(-(X-m3).^2/2/c5))/s5;

if ~exist('syms')
    warning('I can''t draw the remaining figure, because you need to have installed the symbolic toolbox')
    return;
end

%We use symbolic calculus to obtain the values associated with the other
%two regions
syms xs ms ts ps cs
ms=m1;
cs=c2;
ps=(m2-m1);
a=(s2-s1)*int(exp(-(xs-(ms+ts*ps))^2/2/cs),'ts',0,1)/sqrt(2*pi*cs)/s5;
P2=subs(a,'xs',X);

syms xs ms ts ps cs
ms=m2;
cs=c4;
ps=m3-m2;
a=(s4-s3)*int(exp(-(xs-(ms+ts*ps))^2/2/cs),'ts',0,1)/sqrt(2*pi*cs)/s5;
P4=subs(a,'xs',X);
 
P=P1+P2+P3+P4+P5;
figure
plot(X,P,'LineWidth',2,'Color','g');
title('Densidad de Probabilidad','FontSize',12);
xlabel('Concentracion (%) ','FontSize',12);
ylabel('p(Concentracion)','FontSize',12);
hold on;


