%% Robust Regression Using Gibbs Sampler
% Consider a Bayesian linear regression model containing one predictor, and a t distributed disturbance variance with a profiled degrees of freedom parameter \nu.
% \lambda_j\sim IG(\nu/2,2/\nu).
% \varepsilon_j\vert\lambda_j\sim N(0,\lambda_j\sigma^2)
% f(\beta,\sigma^2)\propto \frac{1}{\sigma^2}
% These assumptions imply:
% \varepsilon_j\sim t(0,\sigma^2,\nu)
% \lambda_j\vert\varepsilon_j\sim IG\left(\frac{\nu + 1}{2},\frac{2}{\nu + \varepsilon_j^2/\sigma^2}\right)
% \lambda is a vector of latent scale parameters that attributes low precision to observations far from the regression line.  \nu is a hyperparameter controlling the influence of \lambda on the observations.
% For this problem, the Gibbs sampler is well suited to estimate the coefficients because you can simulate the parameters of a Bayesian linear regression model conditioned on \lambda, and then simulate \lambda from its conditional distribution.
% Generate n = 100 responses from y_t = 1 + 2x_t + e_t, where x \in [0,2] and e_t\sim N(0,0.5^2).

rng('default');
n = 100;
x = linspace(0,2,n)';
b0 = 1;
b1 = 2;
sigma = 0.5;
e = randn(n,1);
y = b0 + b1*x + sigma*e;

%
% Introduce outlying responses by inflating all responses below x = 0.25 by a factor of 3.

y(x < 0.25) = y(x < 0.25)*3;

%
% Fit a linear model to the data. Plot the data and the fitted regression line.

Mdl = fitlm(x,y)
figure;
plot(Mdl);
hl = legend;
hold on;

% The simulated outliers appear to influence the fitted regression line.
%
% Implement this Gibbs sampler:
% Draw parameters from the posterior distribution of \beta,\sigma^2\vert y,x,\lambda. Deflate the observations by \lambda, create a diffuse prior model with two regression coefficients, and draw a set of parameters from the posterior. The first regression coefficient corresponds to the intercept, so specify that bayeslm not include an intercept.
% Compute residuals.
% Draw values from the conditional posterior of \lambda.
% Run the Gibbs sampler for 20,000 iterations and apply a burn-in period of 5,000. Specify \nu = 1, preallocate for the posterior draws, and initialize \lambda to a vector of ones.

m = 20000;
nu = 1;
burnin = 5000;
lambda = ones(n,m + 1);
estBeta = zeros(2,m + 1);
estSigma2 = zeros(1,m + 1);
for j = 1:m
    yDef = y./sqrt(lambda(:,j));
    xDef = [ones(n,1) x]./sqrt(lambda(:,j));
    PriorMdl = bayeslm(2,'Model','diffuse','Intercept',false);
    [estBeta(:,j + 1),estSigma2(1,j + 1)] = simulate(PriorMdl,xDef,yDef);
    ep = y - [ones(n,1) x]*estBeta(:,j + 1);
    sp = (nu + 1)/2;
    sc =  2./(nu + ep.^2/estSigma2(1,j + 1));
    lambda(:,j + 1) = 1./gamrnd(sp,sc);
end

% A good practice is to diagnose the MCMC sampler by examining trace plots. For brevity, this example skips this task.
%
% Compute the mean of the draws from the posterior of the regression coefficients. Remove the burn-in period draws.

postEstBeta = mean(estBeta(:,(burnin + 1):end),2)

% The estimate of the intercept is lower and the slope is higher than the estimates returned by fitlm.
%
% Plot the robust regression line with the regression line fitted by least squares.

h = gca;
xlim = h.XLim';
plotY = [ones(2,1) xlim]*postEstBeta;
plot(xlim,plotY,'LineWidth',2);
hl.String{4} = 'Robust Bayes';

% The regression line fit using robust Bayesian regression appears to be a better fit.
% Copyright 2012 The MathWorks, Inc.