function S = demo_mcmc(nr,sg)

 
if nargin < 2
   sg=4;
end
if nargin < 1
   nr =1;
end

N=10000;
if nr == 2
mu1=-4;
mu2=4;
sg1=2;
sg2=2;
p1=.5;
p2=.5;
x = -12:0.1:12;
y1 = p1*Expo(x,mu1,sg1)+p2*Expo(x,mu2,sg2);
y1_max=max(y1);
y1=y1./y1_max;
ax=[-12 12 0 1];
RN = rand(1,N);
I1=find(RN<=p1);
I2=find(RN>p1);
LI1 = length(I1);
LI2 = length(I2);
X1 = sg1*randn(1,LI1)+mu1;
X2 = sg2*randn(1,LI2)+mu2;
X(I1)=X1;
X(I2)=X2;
elseif nr == 1
mu1=0;
sg1=1;
x = -5:0.1:5;
y1 = Expo(x,mu1,sg1);
y1_max=max(y1);
y1=y1./y1_max; 
ax=[-5 5 0 1];
X = sg1*randn(1,N) + mu1;
end

accept=0;
S=zeros(1,N);
for i=1:N
 i  
 figure(1);
 subplot(2,1,1);
   plot(x,y1,'y','Linewidth',2);hold on;
   y2 = Expo(x,S(i),sg);
   y2_max=max(y2);
   y2=y2./y2_max;
   axis(ax);
   plot(x,y2,'r','Linewidth',2);
   axis(ax);
 subplot(2,1,2);  
   plot(x,y1,'y','Linewidth',2);hold on;
   y2 = Expo(x,S(i),sg);
   y2_max=max(y2);
   y2=y2./y2_max;
   axis(ax);

   Ts = S(i) + sg*randn;
   if nr == 2
   F = (p1*Expo(Ts,mu1,sg1)+p2*Expo(Ts,mu2,sg2))/...
      (p1*Expo(S(i),mu1,sg1)+p2*Expo(S(i),mu2,sg2));
   elseif nr == 1
   F = (Expo(Ts,mu1,sg1))/...
       (Expo(S(i),mu1,sg1));
   end
   A = min(1,F);
   r = rand;
   if r <= A
      S(i+1)=Ts;
      accept=accept+1;
   else
      S(i+1)=S(i);
   end
  
   [M,C]=hist(S(2:i+1),30); 
   [Mx,Cx]=hist(X(1:i),30);
   M_max=max(M);
   subplot(2,1,1);
   bar(C,M./M_max,1);
   Mx_max=max(Mx);
   subplot(2,1,2);
   bar(Cx,Mx./Mx_max,1);
   
   subplot(2,1,1);
   if nr == 2
      plot([-12,-12],[0,accept/i],'g','linewidth',10);
   elseif nr == 1
      plot([-5,-5],[0,accept/i],'g','linewidth',10);

   end
   subplot(2,1,1);
   drawnow;hold off;
   subplot(2,1,2);
   drawnow; hold off;
   
   
end