function S = demo_Gibbs(sg1,sg2,theta)

% function S = demo_Gibbs(sg1,sg2,theta)

if nargin < 3
   theta = pi/4;
end
if nargin < 1
   sg1=10;
   sg2=0.6;
end
   N = 10000;

Sig = [sg1^2 0 ; 0 sg2^2];
R = [cos(theta) -sin(theta) ; sin(theta) cos(theta)];
Sig = R*Sig*R';
sSig1 = sqrt(Sig(1,1)-Sig(1,2)^2/Sig(2,2));
sSig2 = sqrt(Sig(2,2)-Sig(1,2)^2/Sig(1,1));
[U,D,V]=svd(Sig);
K = U*sqrt(D);
X = K * randn(2,N);

S=zeros(2,N);
for i=1:N
   
   figure(1);
   subplot(2,1,1);
   plotGauss(0,0,Sig(1,1),Sig(2,2),Sig(1,2),'r');hold on;
   axis equal;
   subplot(2,1,2);
   plotGauss(0,0,Sig(1,1),Sig(2,2),Sig(1,2),'r');hold on;
   axis equal;

   mu1 = (Sig(1,2)/Sig(2,2))*S(2,i);
   S(1,i+1)=sSig1*randn + mu1;
   mu2 = (Sig(1,2)/Sig(1,1))*S(1,i+1);
   S(2,i+1)=sSig2*randn + mu2;
   
   subplot(2,1,1);
   plot(S(1,2:i+1),S(2,2:i+1),'y+','linewidth',2);hold off;
   subplot(2,1,2);
   plot(X(1,1:i),X(2,1:i),'g+','linewidth',2);hold off;

   drawnow;
end