% T-61.5140 Machine Learning: Advanced Probabilistic Methods
% Solutions for the 8th exercise session (4th and 11th of April, 2008)
% Tapani Raiko

% Problem 1(a)
figure(1); clf;
xx=-5:0.01:5;
k=1.5;
yy=k*gaussian(xx,0,2);
for i=1:2,
  nsamples = 100*(i==1) + 1000*(i==2);
  x = randn(1,nsamples)*sqrt(2);
  y = k*rand(1,nsamples).*gaussian(x,0,2);
  accepted = find(y<unknown_p(x));
  subplot(1,2,i);
  plot(xx,[unknown_p(xx);yy]); 
  hold on;
  plot(x,y,'.c'); 
  plot(x(accepted),y(accepted),'.k');
end;

% Problem 2(a)
figure(2); clf;
importance = unknown_p(xx)./gaussian(xx,0,2);
plot(xx,[importance; unknown_p(xx); gaussian(xx,0,2)]);
legend('importance p_u/q','unknown p_u','proposal q');

% accuracy comparison 1(b),2(b),2(c)
% repeat the estimation a thousand times
for i=1:1000,
  nsamples = 100;
  x = randn(1,nsamples)*sqrt(2);
  y = k*rand(1,nsamples).*gaussian(x,0,2);
  accepted = find(y<unknown_p(x));
  importance = unknown_p(x)./gaussian(x,0,2);
  Ex1(i) = mean(x(accepted));
  Etanhx1(i) = mean(tanh(x(accepted)));
  Ex2(i) = sum(x.*importance)/sum(importance);
  Etanhx2(i) = sum(tanh(x).*importance)/sum(importance);
end;
mean(Ex1)
mean(Ex2)
std(Ex1)
std(Ex2)
mean(Etanhx1)
mean(Etanhx2)
std(Etanhx1)
std(Etanhx2)

% Problem 3(b)
figure(3); clf; hold on;
x1(1) = 0;
x2(1) = 0;
for iter=1:50,
  x1(iter+1) = 10/11*x2(iter)+randn(1)*sqrt(1/11);
  plot([x1(iter);x1(iter+1)],[x2(iter);x2(iter)],'b');
  x2(iter+1) = x1(iter+1)+randn(1)*sqrt(0.1);
  plot([x1(iter+1);x1(iter+1)],[x2(iter);x2(iter+1)],'b');
  plot([x1(iter);x1(iter+1)],[x2(iter);x2(iter+1)],'r');
end;

% Problem 4
figure(4); clf;
x3 = 1.0;
for i=1:8,
  % variance of the proposed jump:
  propvar = 0.01*(mod(i,4)==1)+0.1*(mod(i,4)==2)+1.0*(mod(i,4)==3)+10.0*(mod(i,4)==0);
  nsamples = 50*(i<=4)+500*(i>=5);
  subplot(2,4,i); axis([-4 4 -4 4]); hold on;
  title(sprintf('nsamples = %d, prop var = %3.2f',nsamples,propvar))
  x1(1) = 0;
  x2(1) = 0;
  for iter = 1:nsamples,
    prob_old = unknown_p(x1(iter))*unknown_p(x2(iter))*gaussian(x3,x1(iter)+x2(iter),1.0);
    x1prop = x1(iter)+randn(1)*sqrt(propvar);
    x2prop = x2(iter)+randn(1)*sqrt(propvar);
    prob_prop = unknown_p(x1prop)*unknown_p(x2prop)*gaussian(x3,x1prop+x2prop,1.0);
    if prob_prop/prob_old > rand(1), % accept proposed jump?
      x1(iter+1)=x1prop;
      x2(iter+1)=x2prop;
      plot([x1(iter);x1(iter+1)],[x2(iter);x2(iter+1)],'b-');
    else,
      x1(iter+1)=x1(iter);
      x2(iter+1)=x2(iter);
      plot(x1prop,x2prop,'rx');    
    end;
  end;
end;    


