function post = bayesian_anova(X,Y);

% function post = bayesian_anova(X,Y);
%
% X = Design Matrix
% Y = Data

seed = RandStream.create('mt19937ar','seed',sum(100*clock));
RandStream.setDefaultStream(seed);

burnin = 5000;
thinn = 2500;%10000;
nsamp = 250;%200;
nchain = 3;

n = length(Y);
[nx,k] = size(X);

XtX = X'*X;
iXtX = XtX\eye(k);
Bh = XtX\(X'*Y);
Vi = ones(k,1);

nu = n-k;
s2 = (1/(n-k))*(Y-X*Bh)'*(Y-X*Bh);

%num = (nu*s2/2);
%den = (1+nu/2);

burn = zeros(burnin,nchain,k+1);
thin = zeros(thinn,nchain,k+1);
post = zeros(nsamp,nchain,k+1);
Bt = zeros(k,nchain);
st = zeros(nchain);

for c = 1:nchain
    Bt = mvnrnd(Bh',10*s2*iXtX,1)';
    burn(1,c,1:k) = Bt;
    st = sample_sigma(nu,Bt,X,Y);
    burn(1,c,k+1) = st;
end

disp('beginning burn in')

for n = 2:burnin
    for c = 1:nchain
        Vb = burn(n-1,c,k+1)*iXtX;
        Bt = mvnrnd(Bh',Vb)';
        st = sample_sigma(nu,Bt,X,Y);
        burn(n,c,:) = [Bt; st];
    end
end

disp('done with burn in')

post(1,:,:) = burn(end,:,:);
thin(1,:,:) = burn(end,:,:);
%clear burn

for p = 2:nsamp
    for n = 2:thinn
        for c = 1:nchain
            Vb = thin(n-1,c,k+1)*iXtX;
            Bt = mvnrnd(Bh',Vb)';
            st = sample_sigma(nu,Bt,X,Y);
            thin(n,c,:) = [Bt; st];
            if n==thinn
                post(p,c,:) = thin(n,c,:);
                if c==3
                    disp(['step ' num2str(p) ' of ' num2str(nsamp)])
                    thin(1,:,:) = post(p,:,:);
                end
            end
        end
    end
end


function s_samp = sample_sigma(nu,Bh,X,Y);

s2 = (1/nu)*(Y-X*Bh)'*(Y-X*Bh);
m_s = nu*s2/(nu-2);
s_s = sqrt(2*nu^2*s2^2/( (nu-2)^2*(nu-4) ));
lim = [m_s-2.5*s_s m_s+2.5*s_s];
bin = (lim(2)-lim(1))/1000;
s_vec = lim(1):bin:lim(2);
cdf_s = gammainc(s2*nu./(2*s_vec),nu/2,'upper');
r = rand;
[val,idx] = min(abs(r-cdf_s));
s_samp = s_vec(idx);