% hw5b.m
%

load hw5.mat
X = [ones(length(x),1) x];

J = 500;  % number of Gibbs draws
savestep = 1;  % save every (savestep)th draw

[n,k] = size(X);

% matrix to hold gibbs draws
savedraws = zeros(ceil(J/savestep),k);
    
% initialize chain : could choose a different value here
beta = zeros(k,1);
Z = zeros(n,1);

% main Gibbs loop
for j=1:J,

    % save current draw for beta
    saveindex = j/savestep;
	if (saveindex)==round(saveindex),
		savedraws(saveindex,:) = beta';
	end; 

	% draw for latent Z - not a super efficient approach
	for i=1:n,
	    stop = 0;
            while (stop==0),
		Zi = X(i,:)*beta + randn;
		if (Zi > 0) == y(i),
			stop = 1;
		end;
                Z(i) = Zi;
	    end;
	end;

	% draw for beta
    bhat = inv(X'*X)*X'*Z;
	beta = bhat + chol(inv(X'*X))'*randn(k,1);
end;

plot(savedraws)

betadraws=savedraws(101:500,2);
[mean(betadraws) std(betadraws)]
var(betadraws)
hist(betadraws,20)

% end hw5b.m