Documentation

This is machine translation

Translated by Microsoft
Mouseover text to see original. Click the button below to return to the English verison of the page.

Note: This page has been translated by MathWorks. Please click here
To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

conjugateblm

Bayesian linear regression model with conjugate prior for data likelihood

Description

The Bayesian linear regression model object conjugateblm specifies that the joint prior distribution of the regression coefficients and the disturbance variance, that is, (β, σ2) is the dependent, normal-inverse-gamma conjugate model. The conditional prior distribution of β|σ2 is multivariate Gaussian with mean μ and variance σ2V. The prior distribution of σ2 is inverse gamma with shape A and scale B.

The data likelihood is t=1Tϕ(yt;xtβ,σ2), where ϕ(yt;xtβ,σ2) is the Gaussian probability density evaluated at yt with mean xtβ and variance σ2. The specified priors are conjugate for the likelihood, and the resulting marginal and conditional posterior distributions are analytically tractable. For details on the posterior distribution, see Analytically Tractable Posteriors.

In general, when you create a Bayesian linear regression model object, it specifies the joint prior distribution and characteristics of the linear regression model only. That is, the model object is a template intended for further use. Specifically, to incorporate data into the model for posterior distribution analysis, pass the model object and data to the appropriate object function.

Creation

Syntax

PriorMdl = conjugateblm(NumPredictors)
PriorMdl = conjugateblm(NumPredictors,Name,Value)

Description

example

PriorMdl = conjugateblm(NumPredictors) creates a Bayesian linear regression model object (PriorMdl) composed of NumPredictors predictors and an intercept. The joint prior distribution of (β, σ2) is the dependent normal-inverse-gamma conjugate model. PriorMdl is a template defining the prior distributions and dimensionality of β.

example

PriorMdl = conjugateblm(NumPredictors,Name,Value) uses additional options specified by one or more Name,Value pair arguments. Name is a property name, except NumPredictors, and Value is the corresponding value. Name must appear inside single quotes (''). You can specify several Name,Value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Properties

expand all

You can set property values when you create the model object using name-value pair argument syntax, or after model creation using dot notation. For example, to set a more diffuse prior covariance matrix for PriorMdl, a Bayesian linear regression model containing three model coefficients, enter

PriorMdl.V = 100*eye(3);

Number of predictor variables in the Bayesian multiple linear regression model, specified as a nonnegative integer.

NumPredictors must be the same as the number of columns in your predictor data, which you specify during model estimation or simulation.

When specifying NumPredictors, exclude any intercept term for the value.

After creating a model, if you change the value NumPredictors using dot notation, then all these parameters revert to the default values:

  • Variables names (VarNames)

  • The prior mean of β (Mu)

  • The prior covariance matrix of β (V)

Example: 'NumPredictors',2

Data Types: double

Indicate whether to include regression model intercept, specified as the comma-separated pair consisting of 'Intercept' and a value in this table.

ValueDescription
falseExclude an intercept from the regression model. Hence, β is a p-dimensional vector, where p is the value of the NumPredictors property.
trueInclude an intercept in the regression model. Hence, β is a (p + 1)-dimensional vector. MATLAB® prepends a T-by-1 vector of ones to the predictor data during estimation and simulation.

If you include a column of ones in the predictor data for an intercept term, then set Intercept to false.

Example: 'Intercept',false

Predictor variable names for displays, specified as a string vector or cell vector of character vectors. VarNames must contain NumPredictors elements. VarNames(j) is the name of variable in column j of the predictor data set, which you specify during estimation, simulation, and forecasting.

The default is {'Beta(1)','Beta(2),...,Beta(p)}, where p is the value of NumPredictors.

Example: 'VarNames',["UnemploymentRate"; "CPI"]

Data Types: string | cell | char

Mean parameter of the Gaussian prior on β, specified as a numeric scalar or vector.

If Mu is a vector, then it can have NumPredictors or NumPredictors + 1 elements.

  • For NumPredictors elements, conjugateblm sets the prior mean of the NumPredictors predictors only. Predictor variables correspond to the columns in the predictor data set, which is specified during estimation, simulation, or forecasting. conjugateblm ignores the intercept in the model, that is, any intercept retains the default or its current prior mean.

  • For NumPredictors + 1 elements, the first element corresponds to the prior mean of the intercept, and all other elements correspond to the predictor variables.

Example: 'Mu',[1; 0.08; 2]

Data Types: double

Conditional covariance matrix of Gaussian prior on β, specified as a c-by-c symmetric, positive definite matrix. c can be NumPredictors or NumPredictors + 1.

  • If c is NumPredictors, then conjugateblm sets the prior covariance matrix to

    [1e5000V0].

    That is, conjugateblm attributes the default prior covariances to the intercept, and V to the coefficients of the predictor variables in the data. Rows and columns of V correspond to columns (variables) in the predictor data.

  • If c is NumPredictors + 1, then conjugateblm sets the entire prior covariance to V. The first row and column correspond to the intercept. All other rows and columns correspond to columns in the predictor data.

The default value is a flat prior. For an adaptive prior, specify diag(Inf(Intercept + NumPredictors,1)). Adaptive priors indicate zero precision, that is, you want the prior distribution to have as little influence on the posterior distribution as possible.

V is the prior covariance of β up to a factor of σ2.

Example: 'V',diag(Inf(3,1))

Data Types: double

Shape parameter of inverse gamma prior on σ2, specified as a numeric scalar.

A must be at least -(Intercept + NumPredictors)/2.

With B held fixed, as A increases, the inverse gamma distribution becomes taller and more concentrated. This characteristic weighs the prior model of σ2 more heavily than the likelihood during posterior estimation.

For the functional form of the inverse gamma distribution, see Analytically Tractable Posteriors.

Example: 'A',0.1

Data Types: double

Scale parameter of inverse gamma prior on σ2, specified as a positive scalar or Inf.

With A held fixed, as B increases, the inverse gamma distribution becomes taller and more concentrated. This characteristic weighs the prior model of σ2 more heavily than the likelihood during posterior estimation.

Example: 'B',5

Data Types: double

Object Functions

estimateFit parameters of Bayesian linear regression model to data
simulateSimulate regression coefficients and disturbance variance of Bayesian linear regression model
forecastForecast responses of Bayesian linear regression model
plotVisualize prior and posterior densities of Bayesian linear regression model parameters
summarizeDistribution summary statistics of Bayesian linear regression model

Examples

expand all

Consider the multiple linear regression model that predicts U.S. real gross national product (GNPR) using a linear combination of industrial production index (IPI), total employment (E), and real wages (WR).

For all , is a series of independent Gaussian disturbances with a mean of 0 and variance .

Assume that the prior distributions are:

  • . is a 4-by-1 vector of means and is a scaled 4-by-4 positive definite covariance matrix.

  • . and are the shape and scale, respectively, of an inverse gamma distribution.

These assumptions and the data likelihood imply a normal-inverse-gamma conjugate model.

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

p = 3;
Mdl = bayeslm(p,'ModelType','conjugate')
Mdl = 

 
           |  Mean     Std            CI95         Positive       Distribution     
-----------------------------------------------------------------------------------
 Intercept |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(1)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(2)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Beta(3)   |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

Mdl is a conjugateblm Bayesian linear regression model object representing the prior distribution of the regression coefficients and disturbance variance. At the command window, bayeslm displays a summary of the prior distributions.

You can set writable property values of created models using dot notation. Specify the regression coefficient names to the corresponding variable names.

Mdl.VarNames = ["IPI" "E" "WR"]
Mdl = 

 
           |  Mean     Std            CI95         Positive       Distribution     
-----------------------------------------------------------------------------------
 Intercept |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

This example is based on Create Normal-Inverse-Gamma Conjugate Prior Model.

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',["IPI" "E" "WR"]);

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 and .

PosteriorMdl = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood     -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

PosteriorMdl is a conjugateblm model object storing the joint marginal posterior distribution of and given the data. estimate displays a summary of the marginal posterior distributions to the command window. Rows of the summary correspond to regression coefficients and the disturbance variance, and columns to characteristics of the posterior distribution. The characteristics include:

  • CI95, which contains the 95% Bayesian equitailed credible intervals for the parameters. For example, the posterior probability that the regression coefficient of WR is in [1.782, 3.154] is 0.95.

  • Positive, which contains the posterior probability that the parameter is greater than 0. For example, the probability that the intercept is greater than 0 is 0.003.

  • Distribution, which contains descriptions of the posterior distributions of the parameters. For example, the marginal posterior distribution of IPI is t with a mean of 4.39, a standard deviation of 0.14, and 68 degrees of freedom.

Access properties of the posterior distribution using dot notation. For example, display the marginal posterior means by accessing the Mu property.

PosteriorMdl.Mu
ans =

  -24.2494
    4.3913
    0.0011
    2.4683

This example is based on Create Normal-Inverse-Gamma Conjugate Prior Model.

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

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',["IPI" "E" "WR"])
PriorMdl = 

 
           |  Mean     Std            CI95         Positive       Distribution     
-----------------------------------------------------------------------------------
 Intercept |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

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 conditional posterior distributions of given and the data.

[Mdl,condPostMeanBeta,CondPostCovBeta] = estimate(PriorMdl,X,y,...
    'Sigma2',2);
Method: Analytic posterior distributions
Conditional variable: Sigma2 fixed at   2
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive     Distribution    
--------------------------------------------------------------------------------
 Intercept | -24.2494  1.8695  [-27.914, -20.585]    0.000   N (-24.25, 1.87^2) 
 IPI       |   4.3913  0.0301   [ 4.332,  4.450]     1.000   N (4.39, 0.03^2)   
 E         |   0.0011  0.0001   [ 0.001,  0.001]     1.000   N (0.00, 0.00^2)   
 WR        |   2.4683  0.0743   [ 2.323,  2.614]     1.000   N (2.47, 0.07^2)   
 Sigma2    |    2       0       [ 2.000,  2.000]     1.000   Fixed value        
 

estimate returns the 4-by-1 vector of means and the 4-by-4 covariance matrix of the conditional posterior distribution of the in condPostMeanBeta and CondPostCovBeta, respectively. Also, estimate displays a summary of the conditional posterior distribution of .

Display Mdl.

Mdl
Mdl = 

 
           |  Mean     Std            CI95         Positive       Distribution     
-----------------------------------------------------------------------------------
 Intercept |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 IPI       |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 E         |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 WR        |  0      70.7107  [-141.273, 141.273]    0.500   t (0.00, 57.74^2,  6) 
 Sigma2    | 0.5000   0.5000    [ 0.138,  1.616]     1.000   IG(3.00,    1)        
 

Because estimate computes the conditional posterior distribution, it returns the original prior model, not the posterior, in the first position of the output argument list.

This example is based on Estimate Marginal Posterior Distributions.

Create and estimate the marginal posterior distributions. Extract posterior distribution summary statistics for .

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',["IPI" "E" "WR"]);

load Data_NelsonPlosser
X = DataTable{:,PriorMdl.VarNames(2:end)};
y = DataTable{:,'GNPR'};

[PosteriorMdl,estBeta,estBetaCov] = estimate(PriorMdl,X,y);
Method: Analytic posterior distributions
Number of observations: 62
Number of predictors:   4
Log marginal likelihood     -259.348
 
           |   Mean      Std          CI95        Positive       Distribution      
-----------------------------------------------------------------------------------
 Intercept | -24.2494  8.7821  [-41.514, -6.985]    0.003   t (-24.25, 8.65^2, 68) 
 IPI       |   4.3913  0.1414   [ 4.113,  4.669]    1.000   t (4.39, 0.14^2, 68)   
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000   t (0.00, 0.00^2, 68)   
 WR        |   2.4683  0.3490   [ 1.782,  3.154]    1.000   t (2.47, 0.34^2, 68)   
 Sigma2    |  44.1347  7.8020   [31.427, 61.855]    1.000   IG(34.00, 0.00069)     
 

Suppose that if the coefficient of real wages is below 2.5, then a policy is enacted. Although the posterior distribution of WR is known, and so you can calculate probabilities directly, you can estimate the probability using Monte Carlo simulation instead.

Draw 1e6 samples from the marginal posterior distribution of .

NumDraws = 1e6;
rng(1);
BetaSim = simulate(PosteriorMdl,'NumDraws',NumDraws);

BetaSim is a 4-by- 1e6 matrix containing the draws. Rows correspond to the regression coefficient and columns to successive draws.

Isolate the draws corresponding to the coefficient of real wages, and then identify which draws are less than 2.5.

isWR = PosteriorMdl.VarNames == "WR";
wrSim = BetaSim(isWR,:);
isWRLT2p5 = wrSim < 2.5;

Find the marginal posterior probability that the regression coefficient of WR is below 2.5 by computing the proportion of draws that are less than 2.5.

probWRLT2p5 = mean(isWRLT2p5)
probWRLT2p5 =

    0.5362

The posterior probability that the coefficient of real wages is less than 2.5 is about 0.54.

The marginal posterior distribution of the coefficient of WR is a , but centered at 2.47 and scaled by 0.34. Directly compute the posterior probability that the coefficient of WR is less than 2.5.

center = estBeta(isWR);
stdBeta = sqrt(diag(estBetaCov));
scale = stdBeta(isWR);
t = (2.5 - center)/scale;
dof = 68;
directProb = tcdf(t,dof)
directProb =

    0.5361

The posterior probabilities are nearly identical.

This example is based on Estimate Marginal Posterior Distributions.

Create and estimate the marginal posterior distributions. Hold out the last 10 periods of data from estimation so that they can be used for forecasting real GNP. Turn the estimation display off.

p = 3;
PriorMdl = bayeslm(p,'ModelType','conjugate','VarNames',["IPI" "E" "WR"]);

load Data_NelsonPlosser
fhs = 10; % Forecast horizon size
X = DataTable{1:(end - fhs),PriorMdl.VarNames(2:end)};
y = DataTable{1:(end - fhs),'GNPR'};
XF = DataTable{(end - fhs + 1):end,PriorMdl.VarNames(2:end)}; % Future predictor data
yFT = DataTable{(end - fhs + 1):end,'GNPR'};                  % True future responses

PosteriorMdl = estimate(PriorMdl,X,y,'Display',false);

Forecast responses using the posterior predictive distribution and using the future predictor data XF. Plot the true values of the response and the forecasted values.

yF = forecast(PosteriorMdl,XF);

figure;
plot(dates,DataTable.GNPR);
hold on
plot(dates((end - fhs + 1):end),yF)
h = gca;
hp = patch([dates(end - fhs + 1) dates(end) dates(end) dates(end - fhs + 1)],...
    h.YLim([1,1,2,2]),[0.8 0.8 0.8])
hp = 
  Patch with properties:

    FaceColor: [0.8000 0.8000 0.8000]
    FaceAlpha: 1
    EdgeColor: [0 0 0]
    LineStyle: '-'
        Faces: [1 2 3 4]
     Vertices: [4x2 double]

  Show all properties

uistack(hp,'bottom');
legend('True GNPR','Forecasted GNPR','Forecast Horizon','Location','NW')
title('Real Gross National Product: 1909 - 1970');
ylabel('rGNP');
xlabel('Year');
hold off

yF is a 10-by-1 vector of future values of real GNP corresponding to the future predictor data.

Estimate the forecast root mean squared error (RMSE).

frmse = sqrt(mean((yF - yFT).^2))
frmse = 25.5397

Forecast RMSE is a relative measure of forecast accuracy. Specifically, you estimate several models using different assumptions. The model with the lowest forecast RMSE is the best performing model of the ones being compared.

Definitions

expand all

Algorithms

You can reset all model properties using dot notation, for example, PriorMdl.V = diag(Inf(3,1)). For property resets, conjugateblm does minimal error checking of values. Minimizing error checking has the advantage of reducing overhead costs for Markov chain Monte Carlo simulations, which results in efficient execution of the algorithm.

Alternatives

You can also create a Bayesian linear regression model with a conjugate prior using bayeslm.

Introduced in R2017a

Was this topic helpful?