clc
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Gibbs sampling to identify the joint probability of A and B
%p(A=0/B=0)=1/5
%p(A=1/B=0)=4/5

%p(A=0/B=1)=3/5
%p(A=1/B=1)=2/5

%p(B=0/A=0)=1/4
%p(B=1/A=0)=3/4

%p(B=0/A=1)=2/3
%p(B=1/A=1)=1/3
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PAgivenB0=[1/5 1];
PAgivenB1=[3/5 1];
PBgivenA0=[1/4 1];
PBgivenA1=[2/3 1];
vect=[1 1];
COL=[];
COL=[COL; vect];
TI=[50,100,500,1000,5000,10000];
temp1=50;
count=1;
figure
for i=1:1:10000
     B=vect(2);
     if(B==1)
           if(rand<PAgivenB1(1))
            vect(1)=0;
         else
             vect(1)=1;
         end
     else
         temp=PAgivenB0-rand;
         if(rand<PAgivenB0(1))
            vect(1)=0;
         else
             vect(1)=1;
         end
     end
     COL=[COL;vect];

     A=vect(1);
     if(A==1)
          if(rand<PBgivenA1(1))
            vect(2)=0;
         else
             vect(2)=1;
         end
     else
         temp=PBgivenA0-rand;
         if(rand<PBgivenA0(1))
            vect(2)=0;
         else
             vect(2)=1;
         end
     end
     COL=[COL;vect];
     
%      if i==temp1
%          S1=num2str(temp1);
%          S2='Iterations';
%          CS=strcat(S1,{' '},S2);
%          subplot(2,3,count)
%          hist3(COL./length(COL),[2,2])
%          title(CS);
%          xlabel('A');
%          ylabel('B');
%          count=count+1;
%          if count<=6
%              temp1=TI(count);
%          end
%      end
     
     PA0(i)=length(find(COL(:,1)==0))/length(COL);
     PA1(i)=length(find(COL(:,1)==1))/length(COL);
     
     PB0(i)=length(find(COL(:,2)==0))/length(COL);
     PB1(i)=length(find(COL(:,2)==1))/length(COL);
     
     PA0_PB0(i)=length(find(COL(:,1)==0&COL(:,2)==0))/length(COL);
     PA0_PB1(i)=length(find(COL(:,1)==0&COL(:,2)==1))/length(COL);
     PA1_PB0(i)=length(find(COL(:,1)==1&COL(:,2)==0))/length(COL);
     PA1_PB1(i)=length(find(COL(:,1)==1&COL(:,2)==1))/length(COL);
     
end
L1=(length(COL)-length(nonzeros(sum((COL'-repmat([0  0]',1,length(COL))).^2))))/length(COL);
L2=(length(COL)-length(nonzeros(sum((COL'-repmat([1  0]',1,length(COL))).^2))))/length(COL);
L3=(length(COL)-length(nonzeros(sum((COL'-repmat([0  1]',1,length(COL))).^2))))/length(COL);
L4=(length(COL)-length(nonzeros(sum((COL'-repmat([1  1]',1,length(COL))).^2))))/length(COL);

figure
subplot(221)
plot(PA0_PB0)
axis([0 length(PA0_PB0) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=0,B=0)')
subplot(222)
plot(PA0_PB1)
axis([0 length(PA0_PB1) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=0,B=1)')
subplot(223)
plot(PA1_PB0)
axis([0 length(PA1_PB0) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=1,B=0)')
subplot(224)
plot(PA1_PB1)
axis([0 length(PA1_PB1) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=1,B=1)')

figure
subplot(221)
plot(PA0)
axis([0 length(PA0) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=0)')
subplot(222)
plot(PA1)
axis([0 length(PA1) 0 1])
xlabel('Number of Iterations')
ylabel('P(A=1)')
subplot(223)
plot(PB0)
axis([0 length(PB0) 0 1])
xlabel('Number of Iterations')
ylabel('P(B=0)')
subplot(224)
plot(PB1)
axis([0 length(PB1) 0 1])
xlabel('Number of Iterations')
ylabel('P(B=1)')