Replacing Removed Syntaxes of estimate

The estimate function of the Bayesian linear regression models conjugateblm, semiconjugateblm, diffuseblm, empiricalblm, and customblm returns only an estimated model and an estimation summary table. The previous and current syntaxes for estimate are:

[PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary] = estimate(...); % Now issues error 
[PosteriorMdl,Summary] = estimate(...);                                           % Supported

Starting in R2019b, estimate returns the estimated model object in the first output argument position and the estimation summary table in the second output argument position. If you specify output arguments in subsequent positions, then estimate issues this error:

Too many output arguments.

To avoid the error, update your code by following this procedure.

  1. Search your code for instances in which estimate returns any output argument positions after the first.

  2. In the identified instances, determine whether estimate returns marginal or conditional distribution estimates. If the call to estimate uses the 'Beta' or 'Sigma2' name-value pair argument, the function returns conditional distribution estimates. Otherwise, it returns marginal distribution estimates.

  3. For instances returning marginal distribution estimates:

    1. Return PosteriorMdl and Summary when you call estimate.

    2. Estimate the marginal posterior statistics using the appropriate code in this table.

      Output ArgumentPosteriorMdl Is conjugateblm ObjectPosteriorMdl Is empiricalblm Object
      estBeta

      estBeta = PosteriorMdl.Mu;
      or
      estBeta = Summary.Mean(1:(end - 1));

      estBeta = Summary.Mean(1:(end - 1));
      or
      estBeta = mean(PosteriorMdl.BetaDraws,2);

      EstBetaCov

      EstBetaCov = Summary.Covariances(1:(end - 1),1:(end - 1));
      or
      EstBetaCov = PosteriorMdl.V/PosteriorMdl.B/(PosteriorMdl.A - 1);

      EstBetaCov = Summary.Covariances(1:(end - 1),1:(end - 1));
      or 
      EstBetaCov = cov(PosteriorMdl.BetaDraws');

      estSigma2

      estSigma2 = Summary.Mean(end);
      or
      estSigma2 = 1/PosteriorMdl.B/(PosteriorMdl.A - 1);

      estSigma2 = Summary.Mean(end);
      or
      estSigma2 = mean(PosteriorMdl.Sigma2Draws);

      estSigma2Var

      estSigma2Var = Summary.Covariances(end,end);
      or
      estSigma2Var = 1/(PosteriorMdl.B^2*(PosteriorMdl.A - 1)^2*(PosteriorMdl.A - 2));

      estSigma2Var = Summary.Covariances(end,end);
      or
      estSigma2Var = var(PosteriorMdl.Sigma2Draws);

      For examples, see Replace Removed Syntax When Estimating Analytical Marginal Posterior and Replace Removed Syntax When Estimating Numerical Marginal Posterior.

  4. For instances returning conditional distribution estimates, you must return the estimation summary table Summary in the second output argument position. This table describes how to extract conditional posterior estimates from Summary.

    Output ArgumentEstimation
    estBeta

    estBeta = Summary.Mean(1:end – 1);

    EstBetaCov

    EstBetaCov = Summary.Covariances(1:end – 1,1:end – 1)

    estSigma2

    estSigma2 = Summary.Mean(end)

    estSigma2Var

    estSigma2Var = Summary.Covariances(end,end)

    For an example, see Replace Removed Syntax When Estimating Conditional Posterior.

Replace Removed Syntax When Estimating Analytical Marginal Posterior

Open Live Script

This example shows how to replace the removed syntax of estimate when it returns an analytical marginal posterior.

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

  • The intercept and three regression coefficients are random variables with a multivariate normal prior distribution conditional on the disturbance variance. The prior mean is a 4-D vector of zeros, and the prior covariance matrix is the 4-by-4 identity matrix scaled by 10,000.

  • The disturbance variance is a random variable with an inverse gamma prior distribution. The shape and scale parameter values are 3 and 1, respectively.

Create a normal-inverse-gamma conjugate prior model for the linear regression parameters. Set the number of predictors p. Set the regression coefficient names to the corresponding variable names.

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'};

Before R2019b, estimate could return up to six outputs, each summarizing the posterior distribution. The previously supported syntax is:

[PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary] = estimate(PriorMdl,X,y);

For R2019b, estimate supports returning only two outputs: the posterior model PosteriorMdl and the estimation summary table Summary. Estimate the marginal posterior distribution by using the updated syntax. Return the posterior model and estimation summary table.

[PosteriorMdl,Summary] = 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)

Compute posterior estimates by using the new procedure.

estBeta = PosteriorMdl.Mu % Posterior mean of coefficients
estBeta = 4×1

  -24.2494
    4.3913
    0.0011
    2.4683


EstBetaCov = Summary.Covariances(1:(end - 1),1:(end - 1)) % Posterior covariance of coefficients
EstBetaCov = 4×4

   77.1246    0.7713   -0.0024    0.5311
    0.7713    0.0200   -0.0000   -0.0295
   -0.0024   -0.0000    0.0000   -0.0001
    0.5311   -0.0295   -0.0001    0.1218


estSigma2 = Summary.Mean(end) % Posterior mean of disturbance variance
estSigma2 = 
44.1347

estSigma2Var = Summary.Covariances(end,end) % Posterior variance of disturbance variance
estSigma2Var = 
60.8709

Replace Removed Syntax When Estimating Numerical Marginal Posterior

Open Live Script

This example shows how to replace the removed syntax of estimate when it returns a numerical marginal posterior.

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

  • The intercept and three regression coefficients are random variables with a multivariate normal prior distribution conditional on the disturbance variance. The prior mean is a 4-D vector of zeros, and the prior covariance matrix is the 4-by-4 identity matrix scaled by 10,000.

  • The disturbance variance is a random variable with an inverse gamma prior distribution. The shape and scale parameter values are 3 and 1, respectively.

Create a normal-inverse-gamma semiconjugate prior model for the linear regression parameters. Specify the number of predictors p. Set the regression coefficient names to the corresponding variable names. 

p = 3;
PriorMdl = bayeslm(p,'ModelType','semiconjugate','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'};

Before R2019b, estimate could return up to six outputs, each summarizing the posterior distribution. The previously supported syntax is:

[PosteriorMdl,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary] = estimate(PriorMdl,X,y);

For R2019b, estimate supports returning only two outputs: the posterior model PosteriorMdl and the estimation summary table Summary. Estimate the marginal posterior distribution by using the updated syntax. Return the posterior model and estimation summary table.

[PosteriorMdl,Summary] = estimate(PriorMdl,X,y);
Method: Gibbs sampling with 10000 draws
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95        Positive  Distribution 
-------------------------------------------------------------------------
 Intercept | -23.9999  9.0499  [-41.852, -6.163]    0.003     Empirical  
 IPI       |   4.3933  0.1445   [ 4.112,  4.677]    1.000     Empirical  
 E         |   0.0011  0.0003   [ 0.000,  0.002]    1.000     Empirical  
 WR        |   2.4696  0.3571   [ 1.764,  3.169]    1.000     Empirical  
 Sigma2    |  46.9242  8.4732   [33.244, 66.309]    1.000     Empirical

Compute posterior estimates by using the new procedure.

estBeta = Summary.Mean(1:(end - 1)) % Posterior mean of coefficients
estBeta = 4×1

  -23.9999
    4.3933
    0.0011
    2.4696


EstBetaCov = Summary.Covariances(1:(end - 1),1:(end - 1)) % Posterior covariance of coefficients
EstBetaCov = 4×4

   81.9002    0.8161   -0.0025    0.5843
    0.8161    0.0209   -0.0000   -0.0303
   -0.0025   -0.0000    0.0000   -0.0001
    0.5843   -0.0303   -0.0001    0.1275


estSigma2 = Summary.Mean(end) % Posterior mean of disturbance variance
estSigma2 = 
46.9242

estSigma2Var = Summary.Covariances(end,end) % Posterior variance of disturbance variance
estSigma2Var = 
71.7952

Replace Removed Syntax When Estimating Conditional Posterior

Open Live Script

This example shows how to replace the removed syntax of estimate when it returns a conditional posterior.

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

  • The intercept and three regression coefficients are random variables with a multivariate normal prior distribution conditional on the disturbance variance. The prior mean is a 4-D vector of zeros, and the prior covariance matrix is the 4-by-4 identity matrix scaled by 10,000.

  • The disturbance variance is a random variable with an inverse gamma prior distribution. The shape and scale parameter values are 3 and 1, respectively.

Create a normal-inverse-gamma conjugate prior model for the linear regression parameters. Specify the number of predictors p. Set the regression coefficient names to the corresponding variable names.

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'};

Before R2019b, estimate could return up to six outputs; the latter five summarize the conditional posterior distribution. The previously supported syntax is: 

[~,estBeta,EstBetaCov,estSigma2,estSigma2Var,Summary] = estimate(PriorMdl,X,y);

For R2019b, estimate supports returning only two outputs. The conditional posterior estimates are in the estimation summary table in the second output argument position.

Estimate the conditional posterior distribution of the regression coefficients given that the disturbance variance is 10. Return the estimation summary table by using the updated syntax.

[~,Summary] = estimate(PriorMdl,X,y,'Sigma2',10);
Method: Analytic posterior distributions
Conditional variable: Sigma2 fixed at  10
Number of observations: 62
Number of predictors:   4
 
           |   Mean      Std          CI95         Positive     Distribution    
--------------------------------------------------------------------------------
 Intercept | -24.2494  4.1803  [-32.443, -16.056]    0.000   N (-24.25, 4.18^2) 
 IPI       |   4.3913  0.0673   [ 4.259,  4.523]     1.000   N (4.39, 0.07^2)   
 E         |   0.0011  0.0002   [ 0.001,  0.001]     1.000   N (0.00, 0.00^2)   
 WR        |   2.4683  0.1661   [ 2.143,  2.794]     1.000   N (2.47, 0.17^2)   
 Sigma2    |    10      0       [10.000, 10.000]     1.000   Fixed value

Compute posterior estimates by using the new procedure.

estBeta = Summary.Mean(1:end - 1) % Posterior mean of coefficients
estBeta = 4×1

  -24.2494
    4.3913
    0.0011
    2.4683


EstBetaCov = Summary.Covariances(1:end - 1,1:end - 1) % Posterior covariance of coefficients
EstBetaCov = 4×4

   17.4748    0.1748   -0.0005    0.1203
    0.1748    0.0045   -0.0000   -0.0067
   -0.0005   -0.0000    0.0000   -0.0000
    0.1203   -0.0067   -0.0000    0.0276

See Also

Objects

Functions