function [logpmarg,lp,llf,lpost]...
=gs_probit_chib(X,y,k,n,r1,r2,mu0,V0,f0,f1,betamat);
ydraw=y; % starting draws for latent data
R = r1+r2;
betamean=mean(betamat,2);
betachibmat=zeros(r2,1); %will collect the pdf evaluations
% of p(beta|y)
betadraw=betamat(:,end); % use last draw of GS as starting draw
thoucount=1000;
% start main loop
i=1;
for i=1:R
%%%%%%%%%%%%%
% draw latent data (Note: could have avoided this by collecting it in the
% original GS ....)
%%%%%%%%%%%%%
ydraw(f0)= tnormrnd(X(f0,:)*betadraw,1,...
-inf*ones(length(f0),1),zeros(length(f0),1));
%next is only difference to Tobit:
ydraw(f1)= tnormrnd(X(f1,:)*betadraw,1,...
zeros(length(f1),1),inf*ones(length(f1),1));
%%%%%%%%%%%%
% draw beta and evaluate posterior at betamean
%%%%%%%%%%%%
V1=inv(inv(V0)+X'*X);
mu1=V1*(inv(V0)*mu0+X'*ydraw);
betadraw=mvnrnd(mu1,V1)'; % we need this to feed into the draw of ystar
betachibmat(i) = mvnpdf(betamean,mu1,V1);
if i== thoucount
i
thoucount=thoucount+1000;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate and add log-posterior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lpost=log(mean(betachibmat));
'posterior evaluation complete'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate log-prior
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lp=log(mvnpdf(betamean',mu0',V0));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% evaluate the log-likelihood function
% at betamean
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
llf=0;
j=1;%
for j=1:n
yi=y(j);
Xi=X(j,:);
mui=Xi*betamean;
llfi=(1-yi)*log(normcdf(-mui,0,1))+yi*log(normcdf(mui,0,1)); %individual LH component
llf=llf+llfi; % add up over individuals
end
logpmarg=lp+llf-lpost;
'chib routine complete'