function [x_pct,xm_pct,xdm_pct] = x_draws_pct_summary(x_draws,pct)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

n_pct = size(pct,2);
n_rep = size(x_draws,2);
n_c = size(x_draws,1);
% Compute posterior means
xm = mean(x_draws')';
xdm_draws = x_draws-repmat(xm,1,n_rep);
xdm_min = min(min(xdm_draws));
xdm_max = max(max(xdm_draws));
ngrid = 200;
grid = linspace(xdm_min,xdm_max,ngrid)';
cdf_xdm_c = NaN(n_c,ngrid);
for ic = 1:n_c;
    for ig = 1:ngrid;
        cdf_xdm_c(ic,ig) = mean(xdm_draws(ic,:) <= grid(ig));
    end;
end;
cdf_xdm = mean(cdf_xdm_c)';
xdm_pct = NaN(n_pct,1);
for i = 1:n_pct;
    [tmp,ii] = min(abs(cdf_xdm-pct(i)));
    xdm_pct(i) = grid(ii);
end;
xm_pct = pctile(xm,pct);
xvec_draws = reshape(x_draws,n_c*n_rep,1);
x_pct = pctile(xvec_draws,pct);

end