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
sp=1;            % scale for proposal distribution
M=5000;          % number of iterations
tt=zeros(M,2);   % variable for saving samples
tt(1,:)=[t1 t2]; % save starting point
rr=0;            % variable for saving number of rejections

% TÄHÄN OMA KOODISI (alle 15 riviä)
% Metropolis sampling here>>>

% Metropolis sampling here<<<

%Rejection rate
rr/M

% Example figure
clf
subplot(2,2,1)
% show the Metropolis 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',1,'MarkerFaceColor','b')
set(gca,'XLim',[-4 4],'YLim',[-4 4])
% contour at 2 times the standard deviation
q=sqrt(eig([1 r;r 1]))*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 beahves when more samples are sampled
% you may need to change the burnin value
burnin=M/5;
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=min(100,M/2);
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, note that it's accuracy improves with longer chains
acorrtime(tt(M/2:end,:),maxlag)