% Program "mcmc"

clc; close all; clear; tic 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% PROMENLJIVE %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

delta=XXX;                    % Parametar za sirinu uniformne raspodele

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

broj_tacaka=2000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tacke=zeros(2,broj_tacaka);
brojac=0;
tacke(:,1)=[0;0];
stara_tacka=tacke(:,1);
for i=2:broj_tacaka
nova_tacka=stara_tacka+delta*(rand(2,1)-0.5);
stara_p=exp(-(stara_tacka(1)^2)*(stara_tacka(2)^2)-stara_tacka(1)^2-stara_tacka(2)^2);
nova_p=exp(-(nova_tacka(1)^2)*(nova_tacka(2)^2)-nova_tacka(1)^2-nova_tacka(2)^2);
    if rand<(nova_p/stara_p)
        stara_tacka=nova_tacka;
        brojac=brojac+1;
    end
tacke(:,i)=stara_tacka;
end

uzorak=tacke(1,:).^2;
integral=sum(uzorak)/broj_tacaka;

num_integral_1=integral2(@(x, y) x.*x.*exp(-x.*x-y.*y-x.*x.*y.*y), -20, 20, -20, 20);
num_integral_2=integral2(@(x, y) exp(-x.*x-y.*y-x.*x.*y.*y), -20, 20, -20, 20);
num_integral=num_integral_1/num_integral_2;

fprintf('\n')
fprintf('-----------------------------------------------------\n')
fprintf('\n')

fprintf ( 1, 'Vrednost integrala dobijena numerickom integracijom na mrezi je:   %f \n', num_integral);
fprintf ( 1, 'Vrednost integrala dobijena Metropolisovom MC metodom je:          %f \n', integral);
fprintf ( 1, 'Procenat prihvacenih tacaka:                                       %2.1f%% \n', brojac/(broj_tacaka-1)*100);

figure
plot(tacke(1,:),tacke(2,:),'-o', ...
    'LineWidth',1, ...
    'MarkerSize',2, ...
    'MarkerEdgeColor','b',...
    'MarkerFaceColor','b')
hold on
axis square
xlabel('\bf \fontsize{15} \fontname{Cambria} x koordinata');
ylabel('\bf \fontsize{15} \fontname{Cambria} y koordinata');
x = linspace(-2,2);
y = linspace(-2,2);
[X,Y] = meshgrid(x,y);
Z = exp(-(X.^2.*Y.^2+X.^2+Y.^2));
contour(X,Y,Z,7,'-k', ...
    'LineWidth',1.2)

fprintf('\n')
fprintf('-----------------------------------------------------\n')
fprintf('\n')

toc