%% Estimate Maximum A Posteriori Probability Using Monte Carlo
% The maximum a posteriori probability (MAP) estimate is the posterior mode, that is, the parameter value that yields the maximum of the posterior pdf. If the posterior is analytically intractable, then you can use Monte Carlo sampling to estimate the MAP.
% Consider the linear regression model in docid:econ_ug.bvj_ani-2.
% Load the Nelson-Plosser data set. Create variables for the response and predictor series.

load Data_NelsonPlosser
varNames = {'IPI' 'E' 'WR'};
X = DataTable{:,varNames};
y = DataTable{:,'GNPR'};

%
% Create a normal-inverse-gamma conjugate prior model for the linear regression parameters. Specify the number of predictors p and the variable names.

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',varNames)

%
% Estimate the marginal posterior distributions of \beta and \sigma^2.

rng(1); % For reproducibility
PosteriorMdl = estimate(PriorMdl,X,y);

% The display includes the marginal posterior distribution statistics.
% Extract the posterior mean of \beta from the posterior model, and extract the posterior covariance of \beta from the estimation summary returned by summarize.

estBetaMean = PosteriorMdl.Mu;
Summary = summarize(PosteriorMdl);
EstBetaCov = Summary.Covariances{1:(end - 1),1:(end - 1)};

% estBetaMean is a 4-by-1 vector representing the mean of the marginal posterior of \beta.  EstBetaCov is a 4-by-4 matrix representing the covariance matrix of the posterior of \beta.
%
% Draw 10,000 parameter values from the posterior distribution.

rng(1); % For reproducibility
[BetaSim,sigma2Sim] = simulate(PosteriorMdl,'NumDraws',1e5);

% BetaSim is a 4-by-10,000 matrix of randomly drawn regression coefficients. sigma2Sim is a 1-by-10,000 vector of randomly drawn disturbance variances.
%
% Transpose and standardize the matrix of regression coefficients. Compute the correlation matrix of the regression coefficients.

estBetaStd = sqrt(diag(EstBetaCov)');
BetaSim = BetaSim';
BetaSimStd = (BetaSim - estBetaMean')./estBetaStd;
BetaCorr = corrcov(EstBetaCov);
BetaCorr = (BetaCorr + BetaCorr')/2; % Enforce symmetry

%
% Because the marginal posterior distributions are known, evaluate the posterior pdf at all simulated values.

betaPDF = mvtpdf(BetaSimStd,BetaCorr,68);
a = 34;
b = 0.00069;
igPDF = @(x,ap,bp)1./(gamma(ap).*bp.^ap).*x.^(-ap-1).*exp(-1./(x.*bp));...
    % Inverse gamma pdf
sigma2PDF = igPDF(sigma2Sim,a,b);

%
% Find the simulated values that maximize the respective pdfs, that is, the posterior modes.

[~,idxMAPBeta] = max(betaPDF);
[~,idxMAPSigma2] = max(sigma2PDF);
betaMAP = BetaSim(idxMAPBeta,:);
sigma2MAP = sigma2Sim(idxMAPSigma2);

% betaMAP and sigma2MAP are the MAP estimates.
%
% Because the posterior of \beta is symmetric and unimodal, the posterior mean and MAP should be the same. Compare the MAP estimate of \beta with its posterior mean.

table(betaMAP',PosteriorMdl.Mu,'VariableNames',{'MAP','Mean'},...
    'RowNames',PriorMdl.VarNames)

% The estimates are fairly close to one another.
%
% Estimate the analytical mode of the posterior of \sigma^2. Compare it to the estimated MAP of \sigma^2.

igMode = 1/(b*(a+1))
sigma2MAP

% These estimates are also fairly close.
%
% Copyright 2018 The MathWorks, Inc.