%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Gibbs sampling example
%    Implements a Gibbs
%     sampler for the 
%     mean and variance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [mu,phi] = GibbsSecondExpl(x, trueMu, trueSigma)

trace=1;


thin = 1;
step = .1;
iterations = 100000;
burnin = 2;

muInit=50;
phiInit =.25;
particlesInit = muInit;

%priorMu = 1;
%priorVar = 1/sigma^2; 

%%Create space for the stored variables
mu = zeros(1,iterations);
phi = zeros(1,iterations);

mu(1) = muInit;
phi(1) = phiInit;

figure(1)
clf
plot(trueMu, 1/trueSigma^2,'r.')
xlabel('\mu','fontsize',14)
ylabel('\phi','fontsize',14)
pause
hold on


tic
for(i=2:iterations)

    %%update the mean (Sample from Normal(\bar x),1/(N\phi))
    proposalMu = mu(i-1) + randn*step;

    %log scale
    proposed    = sum( - phi(i-1)*(x-proposalMu).^2/2);
    currentlyAt = sum( - phi(i-1)*(x-mu(i-1)).^2/2);  


    ratio = proposed-currentlyAt;

    if(log(rand)<ratio)
        mu(i) = proposalMu;
    else
        mu(i) = mu(i-1);
    end



    %%update the Precision (Sample from Gamma(N/2, 0.5*\sum(x-\mu)))
    proposalPhi=-1;
    while(proposalPhi<0)
        proposalPhi = phi(i-1) + randn*step;
    end
    proposed    = -log(proposalPhi)+sum(log( sqrt(proposalPhi)/sqrt(2*pi))  - proposalPhi*(x-mu(i)).^2/2);
    currentlyAt = -log(phi(i-1))+sum(log(sqrt(phi(i-1))/sqrt(2*pi))  - phi(i-1)*(x-mu(i)).^2/2);  



    ratio = proposed-currentlyAt;

    if(log(rand)<ratio)
        phi(i) = proposalPhi;
    else
        phi(i) = phi(i-1);
    end


   if(mod(i,10000)==0)
        figure(1)
        clf
        xlabel('\mu','fontsize',14)
        ylabel('\phi','fontsize',14)
        toc
        i
        hold on
        plot(mu(1:i),phi(1:i));
        plot(trueMu, 1/trueSigma^2,'r.')
        pause(.01);
        tic
   end

end
toc


if(trace==1)
    figure(2)
    clf;
    subplot(2,1,1)
    mu = mu(1:thin:length(mu));         %% thin samples
    plot(mu)

    mu = mu((burnin+1):length(mu));     %% chop off burnin 

    subplot(2,1,2)
    %normalizedHist(mu,.2)
    hist(mu)
    hold on;
    plot([trueMu,trueMu],[0,length(mu)/10],'r','linewidth',2)

    figure(3)
    clf;
    subplot(2,1,1)
    phi = phi(1:thin:length(phi));      %% thin samples
    plot(phi)

    phi = phi((burnin+1):length(phi));  %% chop off burnin

    
    subplot(2,1,2)
    %normalizedHist(phi,.1)
    hist(phi)
    hold on;
    plot([1/trueSigma^2,1/trueSigma^2],[0,length(mu)/10],'r','linewidth',2)
    1/trueSigma^2
end    