%% Tune Slice Sampler For Posterior Estimation
% This example shows how to improve the slice sampler for posterior
% estimation and inference of a Bayesian linear regression model.
%%
% When estimating the posterior composed of the data likelihood and
% semiconjugate or custom prior models, |estimate| uses an MCMC sampler. If
% a trace plot of the sample shows transient behavior or very high serial
% correlation, or you want to store few samples from the posterior, then
% you can specify a burn-in sample or thinning.
%%
% Consider the multiple linear regression model that predicts the US real
% gross national product (|GNPR|) using a linear combination of industrial
% production index (|IPI|), total employment (|E|), and real wages (|WR|).
%
% $$\texttt{GNPR}_t=\beta_0+\beta_1\texttt{IPI}_t+\beta_2\texttt{E}_t+\beta_3\texttt{WR}_t+\varepsilon_t.$$
%
% For all $t$ time points, $\varepsilon_t$ is a series of independent Gaussian
% disturbances with a mean of 0 and variance $\sigma^2$. Assume that the
% prior distributions are:
%
% * $\beta_j\vert\sigma^2$ is 4-dimensional _t_ distribution with 50
% degrees of freedom for each component and the identity matrix for the
% correlation matrix.  Also, the distribution is centered at ${\left[
% {\begin{array}{*{20}{c}} { - 25}&4&0&3 \end{array}} \right]^\prime }$ and
% each component is scaled by the corresponding elements of the vector
% ${\left[ {\begin{array}{*{20}{c}} {1}&1&1&1 \end{array}}
% \right]^\prime }$. 
%
% * $\sigma^2 \sim IG(3,1)$. 
%
% |bayeslm| treats these assumptions and the data likelihood as if the
% corresponding posterior is analytically intractable. 
%%
% Declare a MATLAB(R) function that:
%
% * Accepts values of $\beta$ and $\sigma^2$ together in a column vector,
% and values of the hyperparameters.
% * Returns the value of the joint prior distribution,
% $\pi\left(\beta,\sigma^2\right)$, given the values of $\beta$ and
% $\sigma^2$.
%
% <include>priorMVTIG.m</include>
%
%%
% Create an anonymous function that operates like |priorMVTIG|, but accepts
% only the parameter values, and holds the hyperparameter values fixed
% at their values.
dof = 50;
C = eye(4);
ct = zeros(4,1);
st = ones(4,1);
a = 3;
b = 1;
logPDF = @(params)priorMVTIG(params,ct,st,dof,C,a,b);
%%
% Create a custom joint prior model for the linear regression parameters.
% Specify the number of predictors, |p|. Also, specify the function handle
% for |priorMVTIG| and the variable names.
p = 3;
PriorMdl = bayeslm(p,'ModelType','custom','LogPDF',logPDF,...
    'VarNames',["IPI" "E" "WR"]);
%%
% |PriorMdl| is a |customblm| Bayesian linear regression model object
% representing the prior distribution of the regression coefficients and
% disturbance variance.
%%
% Load the Nelson-Plosser data set.  Create variables for the response and
% predictor series.
load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};  
%%
% Estimate the marginal posterior distributions of $\beta$ and $\sigma^2$.
% To determine an appropriate burn-in period, specify no burn-in period.
% Specify a width for the slice sampler that is close to the posterior
% standard deviation of the parameters assuming a diffuse prior model.
% Draw |15000| samples.
width = [20,0.5,0.01,1,20];
numDraws = 15e3;
rng(1) % For reproducibility
PosteriorMdl1 = estimate(PriorMdl,X,y,'BurnIn',0,'Width',width,...
    'NumDraws',numDraws);
%%
% |PosteriorMdl1| is an |empiricalblm| model object.  The |BetaDraws| and
% |Sigma2Draws| properties store the samples from the slice sampler.
%%
% Draw trace and autocorrelation plots for the
% Markov chains.
figure;
for j = 1:4
    subplot(2,3,j)
    plot(PosteriorMdl1.BetaDraws(j,:));
    axis tight
    xlabel('Sample index')
    title(sprintf(['Trace Plot ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
plot(PosteriorMdl1.Sigma2Draws);
axis tight
xlabel('Sample index')
title(['Trace Plot ' char(8212) ' Sigma2']);

figure;
for j = 1:4
    subplot(2,3,j)
    autocorr(PosteriorMdl1.BetaDraws(j,:));
    axis tight
    title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
autocorr(PosteriorMdl1.Sigma2Draws);
axis tight
title(['Correlogram ' char(8212) ' Sigma2']);
%%
% All parameters, except the intercept, show significant transient effects
% and high autocorrelation.  |estimate| checks for high correlation, and,
% if present, issues warnings.  Also, the parameters do not seem to have
% converged to their stationary distribution.
%%
% Estimate the posterior distribution again.  To reduce transient effects,
% use the default burn-in period (|5000|, which is the default).  To reduce
% the affects of serial correlation, thin by a factor of 20 and specify
% |15e4/thin| effective number of draws.  Also, numerically stabilize the
% estimation of $\sigma^2$ by specifying reparameterize to
% $\log(\sigma^2)$.  Plot trace and autocorrelation plots.
thin = 20;
numDraws = 15e4/thin;
PosteriorMdl2 = estimate(PriorMdl,X,y,'Thin',thin,'Width',width,...
    'Reparameterize',true,'NumDraws',numDraws);

figure;
for j = 1:4
    subplot(2,3,j)
    plot(PosteriorMdl2.BetaDraws(j,:));
    axis tight
    xlabel('Sample index')
    title(sprintf(['Trace Plot ' char(8212) '%s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
plot(PosteriorMdl2.Sigma2Draws);
axis tight
xlabel('Sample index')
title(['Trace Plot ' char(8212) ' Sigma2']);

figure;
for j = 1:4
    subplot(2,3,j)
    autocorr(PosteriorMdl2.BetaDraws(j,:));
    title(sprintf(['Correlogram ' char(8212) ' %s'],PriorMdl.VarNames{j}));
end
subplot(2,3,5)
autocorr(PosteriorMdl2.Sigma2Draws);
title(['Correlogram  ' char(8212) ' Sigma2']);
%%
% The trace plots suggest that the Markov chains:
%
% * Have reached their stationary distributions because the distribution of
% the samples have a relatively constant mean and variance.  Specifically, the
% distribution of the points does not seem to change.
% * Are mixing fairly quickly because the serial correlation is moderate,
% for some parameters, to very low.
% * Do not show any transient behavior.
%
%%
% These results suggest that you can proceed with a sensitivity analysis,
% or with posterior inference.

%% 
% Copyright 2012 The MathWorks, Inc.