function  [theta_draws,resid_draws]=draw_stacked_regression_posterior(Y,X,V,theta_bar, V_theta,ndraws)

% Draw from the posterior of coefficients for a regression using
% informative prior, following the approach of Gelman, Carlin, Stern,
% Dunson, Vehtari, and Rubin (2013)
%
% model is:
% Y= X * theta  + U,   U  is Multi Normal (zeros(T_Star,1), V)
% prior for theta is Multi Normal (theta_bar, V_theta)
%
%  input:
%  Y is T_star x1,  X is T_star x m, V is T_star x T_star
%  theta_bar is mx1,   V_theta is mxm
% output:  theta_draws is mx1
%        or m x ndraws if there is a sixth argument for the number of draws

% Copyright (c) 2013-2018, Francesco Corsello, All Rights Reserved.

b_plurimae=not(nargin<6);

[Tst,m]=size(X);


Y_diam=[Y; theta_bar];   X_diam=[X; eye(m)];


Tm=Tst+m; 


V_diam=sparse(Tm,Tm);
V_diam(1:Tst,1:Tst)=V;  V_diam(Tst+1:end,Tst+1:end)=V_theta;


if isdiag(V_diam)
    
    invch_V_diam=diag(diag(V_diam).^(-0.5));
    X_dd=invch_V_diam*X_diam;  
    Y_dd=invch_V_diam*Y_diam;
    
else
    
    try
        ch_V_diam=chol(V_diam,'lower');
    catch
        ch_V_diam=chol( correct_precise_eig(make_sim_precise( full(V_diam))),'lower');
    end
    
    X_dd=ch_V_diam\X_diam;   
    Y_dd=ch_V_diam\Y_diam;
    
    
end

X_dd_pr_X_dd=X_dd'*X_dd;
X_dd_pr_Y=X_dd'* Y_dd;

var_theta_post=X_dd_pr_X_dd\eye(m);
mean_theta_post=var_theta_post*X_dd_pr_Y;

if b_plurimae
    theta_draws=zeros(m,ndraws);
    parfor j=1:ndraws
        theta_draws(:,j)=transpose(sample_multinorm(  mean_theta_post' , ...
            var_theta_post  ));
    end
    resid_draws=  repmat(Y,[1,ndraws])  -  X * theta_draws    ;
    
else
    theta_draws=transpose(sample_multinorm(  mean_theta_post' , ...
        var_theta_post  ));
    resid_draws=Y- X * theta_draws;
    
end



