function samples = PostSampling(data, dist_type, Ndraws)

switch dist_type
    case 'tmu'
        samples = tMuSampling(data, Ndraws); %posterior for mu
    case 'tpred'
        samples = tPredSampling(data, Ndraws); %posterior predictive
    case 'IG'
        samples = IGSampling(data, Ndraws); %posterior for sigma
end

end

function samples = tMuSampling(data, Ndraws)

NN = length(data); %number of data points

%mean and mean square of log data
Xbar = sum(log(data))./NN;
X2bar = sum(log(data).^2)./NN;

%paramters of t-distribution
mumu = Xbar;
numu = NN-1;
sigmu = sqrt((X2bar - Xbar.^2)./(NN-1));

%calculate pdf values from t-dist
pd = makedist('tLocationScale','mu',mumu,'sigma',sigmu,'nu',numu);

%sample from pdf
samples = random(pd,[Ndraws,1]);

end

function samples = tPredSampling(data, Ndraws)

NN = length(data); %number of data points

%mean and mean square of log data
Xbar = sum(log(data))./NN;
X2bar = sum(log(data).^2)./NN;

%paramters of t-distribution
mutilde = Xbar;
nutilde = NN-1;
sigtilde = sqrt(((NN+1)./(NN-1)).*(X2bar - Xbar.^2));

%calculate pdf values from t-dist
pd = makedist('tLocationScale','mu',mutilde,'sigma',sigtilde,'nu',nutilde);

%sample from pdf
samples = random(pd,[Ndraws,1]);

end

function samples = IGSampling(data, Ndraws)

NN = length(data); %number of data points

%mean and mean square of log data
Xbar = sum(log(data))./NN;
X2bar = sum(log(data).^2)./NN;

%parameters for gamma distribution
asig = (NN-1)./2;
bsig = NN.*(X2bar - Xbar.^2)./2;

%calculate pdf values from gamma distribution
pd = makedist('Gamma','a',asig,'b',1./bsig);

%sample from pdf
samples = random(pd,[Ndraws,1]);
samples = (1./sqrt(samples)); %to change from gamma for precision to inverse gamma for standard deviation

end