clear all
close all
clc

%% Definice parametru a nastaveni Gibbsova vzorkovace
mu = [0
      0];
rho = 0.5; % lze menit pro zkoumani vlastnosti ruzne perzistentniho Gibbsova vzorkovace

%pocet generovanych vzorku + jedna pocatecni podminka
S = 100000+1;
S1 = 90000; %pocet ponechanych vzorku
S0 = S-1-S1; %pocet vyhozenych vzorku

theta = zeros(2,S); %matice pro ukladani vzorku
%pocatecni podminka (lze experimentovat s ruznym nastavenim)
theta(:,1)=[10
            -10];

for s=2:S
   %generujeme theta1|theta2
     mu12 = mu(1)+rho*(theta(2,s-1)-mu(2));
     Sigma12 = 1-rho^2;
    theta(1,s)=mu12+randn*sqrt(Sigma12);
    
   %generujeme theta2|theta1
     mu21 = mu(2)+rho*(theta(1,s)-mu(1));
     Sigma21 = 1-rho^2;
    theta(2,s)= mu21+randn*sqrt(Sigma21);    
end

%% Zobrazeni prvnich 20 vzorku
 figure
 plot(theta(1,1:20),theta(2,1:20),'*-')
 xlabel('\theta_1')
 ylabel('\theta_2')

%odstranime prvnich S0+1 vzorku
theta(:,1:S0-1)=[];

%% Zobrazeni prvnich 20 vzorku (po odstraneni vlivu pocatecnich podminek)
 figure
 plot(theta(1,1:20),theta(2,1:20),'*-')
 xlabel('\theta_1')
 ylabel('\theta_2')
 
%% Vypocet strednich hodnot a rozptylu
% posteriorni momenty
 E_theta = mean(theta,2);
 D_theta = var(theta,0,2);
 
 disp('E      D')
 disp([E_theta D_theta])
 
%% Graficky pohled na generovane vzorky
k = 100; %vezmeme kazdy 100. vzorek
figure
 subplot(2,1,1)
 plot(theta(1,1:k:end))
 ylabel('\theta_1')
 subplot(2,1,2)
 plot(theta(2,1:k:end))
 ylabel('\theta_2')

%% Zobrazeni posteriornich hustot
figure
 subplot(2,1,1)
 hist(theta(1,:),30)
 subplot(2,1,2)
 hist(theta(2,:),30)