[y,Inc,Hsz,WW] = textread('shop.txt','%f %f %f %f');
n=84;
for i=1:n 
    x(i,1)=1; x(i,2)=Inc(i); x(i,3)=Hsz(i); x(i,4)=WW(i); end
X = x(:,[1:4]); beta = [0 0 0 0]'; Lo = -10 .* (1-y); Hi =10 .* y;  T=1000; burn=50;
for ch=1:2 for t=1:T
% truncated normal sample between  Lo and Hi   
    Z = rand_nort(X * beta, ones(size(X * beta)), Lo, Hi);
    sigma=inv(X' * X);     betaMLE = inv(X' * X)* X' * Z;
    beta = rand_MVN(1, betaMLE, sigma)';
for j=1:4   betas(t,j,ch)=beta(j); end 
end
end

figure(1)
beta0=betas(burn:T,1,1:2);  subplot(4,2,1), hist(beta0,50)
beta1=betas(burn:T,2,1:2);  subplot(4,2,2), hist(beta1,50)
beta2=betas(burn:T,3,1:2);  subplot(4,2,3), hist(beta2,50)
beta3=betas(burn:T,4,1:2);  subplot(4,2,4), hist(beta3,50)
[PSRF] = GRpsrf(betas,T,4,2)
[MPSRF] = GRmpsrf(betas,T,4,2)

