y1=0;y2=0;r=0.8; % parameters of the Normal-distribution
S=[1 r;r 1];     % covariance matrix
t1=-2.5;t2=2.5;  % starting value for the chain
M=1000;           % number of iterations
tt=zeros(M,2);   % variable for saving samples
tt(1,:)=[t1 t2]; % save starting point

% TÄHÄN OMA KOODISI (alle 10 riviä)
% Gibbs sampling here >>>

% Gibbs sampling here <<<

% Example figure
clf
subplot(2,2,1)
% show component-by-component updating of the Gibbs iterations
line(tt(1,1),tt(1,2),'Marker','o')
set(gca,'XLim',[-4 4],'YLim',[-4 4])
line(tt(:,1),tt(:,2))
set(gca,'DataAspectRatio',[1 1 1])

subplot(2,2,2)
% show the iterates from the second half of the sequence
line(tt(M/2:end,1),tt(M/2:end,2),'LineStyle','none','Marker','o','MarkerSize',2,'MarkerFaceColor','b')
set(gca,'XLim',[-4 4],'YLim',[-4 4])
% contour at 2 times the standard deviation
q=sqrt(eig(S))*2;
ellipse(q(2),q(1),pi/4,0,0,'r'); 
set(gca,'DataAspectRatio',[1 1 1])

subplot(2,2,3)
% show how the average behaves when more samples are sampled
burnin=10;
plot(burnin:M,cumsum(tt(burnin:M,:))./repmat([1:length(burnin:M)]',1,2))
line([burnin M],[0 0])
axis tight

subplot(2,2,4)
% show the autocorrelation in the samples in the second half of the sequence
% if you change value of r, you might need to change value of maxlag to get better
% view of maxlag
maxlag=20;
plot(acorr(tt(M/2:end,:),maxlag))
set(gca,'XLim',[1 maxlag],'YLim',[-0.5 1])
line([0 maxlag],[0 0])
% acorrtime can be used to estimate how much you should thin the chain to get
% approximately independent samples
acorrtime(tt(M/2:end,:),maxlag)