# FSRBeda

FSRBeda enables to monitor several quantities in each step of the Bayesian search

## Syntax

• out=FSRBeda(y, X)example
• out=FSRBeda(y, X,Name,Value)example

## Description

 out =FSRBeda(y, X) FSRBeda with all default options.

 out =FSRBeda(y, X, Name, Value) FSRBeda with optional arguments.

## Examples

expand all

### FSRBeda with all default options.

Common part to all examples: load Houses Price Dataset.

load hprice.txt;
% setup parameters
n=size(hprice,1);
y=hprice(:,1);
X=hprice(:,2:5);
[out]=FSRBeda(y,X);

### FSRBeda with optional arguments.

load hprice.txt;
% setup parameters
n=size(hprice,1);
y=hprice(:,1);
X=hprice(:,2:5);
bayes=struct;
n0=5;
bayes.n0=n0;
% set \beta components
beta0=0*ones(5,1);
beta0(2,1)=10;
beta0(3,1)=5000;
beta0(4,1)=10000;
beta0(5,1)=10000;
bayes.beta0=beta0;
% \tau
s02=1/4.0e-8;
tau0=1/s02;
bayes.tau0=tau0;
% R prior settings
R=2.4*eye(5);
R(2,2)=6e-7;
R(3,3)=.15;
R(4,4)=.6;
R(5,5)=.6;
R=inv(R);
bayes.R=R;
[out]=FSRBeda(y,X,'bayes',bayes);

## Related Examples

expand all

### Plot posterior estimates of beta and sigma2.

Plot posterior estimates of beta and sigma2 in the interval (subset size) [20 125].

% In this example for the house price data we monitor the forward plots
% of the parameters of the linear model adding 95% and 99% HPD
% regions. The first 8 panels refer to the elements of $\beta_1$
% and the bottom right-hand panel refer to the estimate of sigma2.
% The horizontal lines correspond to prior values
% The vertical lines refer to the step prior to the introduction of the first outlier
% init = point to start monitoring diagnostics along the FS
% run routine FSRB in order to find the outliers automatically
% setup parameters
n=size(hprice,1);
y=hprice(:,1);
X=hprice(:,2:5);
bayes=struct;
n0=5;
bayes.n0=n0;
% set \beta components
beta0=0*ones(5,1);
beta0(2,1)=10;
beta0(3,1)=5000;
beta0(4,1)=10000;
beta0(5,1)=10000;
bayes.beta0=beta0;
% \tau
s02=1/4.0e-8;
tau0=1/s02;
bayes.tau0=tau0;
% R prior settings
R=2.4*eye(5);
R(2,2)=6e-7;
R(3,3)=.15;
R(4,4)=.6;
R(5,5)=.6;
R=inv(R);
bayes.R=R;
outBA=FSRB(y,X,'bayes',bayes', 'plots',0);
dout=n-length(outBA.ListOut);
% init = initial point to start monitoring
init=20;
xlimL=init; % lower value of xlim
xlimU=125;  % upper value of xlim
outBAeda=FSRBeda(y,X,'bayes',bayes,'init',init, 'conflev', [0.95 0.99]);
% Set font size, line width and line style
figure;
lwd=2.5;
FontSize=14;
linst={'-','--',':','-.','--',':'};
for j=1:5
my_subplot=subplot(3,2,j);
hold('on')
% plot 95% and 99% HPD  trajectories
plot(outBAeda.beta1(:,1),outBAeda.beta1HPD(:,1:2,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','b')
plot(outBAeda.beta1(:,1),outBAeda.beta1HPD(:,3:4,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% plot posterior estimate of beta1_j
plot(outBAeda.beta1(:,1),outBAeda.beta1(:,j+1)','LineStyle',linst{1},'LineWidth',lwd,'Color','k')
% Add the horizontal line which corresponds to prior values
xL = get(my_subplot,'XLim');
line(xL,[beta0(j) beta0(j)],'Color','r','LineWidth',lwd);
% Set ylim
ylimU=max([outBAeda.beta1HPD(:,4,j); beta0(j)]);
ylimL=min([outBAeda.beta1HPD(:,3,j); beta0(j)]);
ylim([ylimL ylimU])
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
% Set xlim
xlim([xlimL xlimU]);
ylabel(['$\hat{\beta_' num2str(j-1) '}$'],'Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
if j>4
xlabel('Subset size m','FontSize',FontSize);
end
end
% Subplot associated with the monitoring of sigma^2
subplot(3,2,6);
%figure()
hold('on')
% 99%
plot(outBAeda.sigma21HPD(:,1),outBAeda.sigma21HPD(:,4:5),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% 95%
plot(outBAeda.sigma21HPD(:,1),outBAeda.sigma21HPD(:,2:3),'LineStyle',linst{2},'LineWidth',lwd,'Color','b')
% Plot 1\/tau1
plot(outBAeda.S21(:,1),1./outBAeda.S21(:,3),'LineWidth',lwd,'Color','k')
ylabel('$\hat{\sigma}^2$','Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
% Set ylim
ylimU=max([outBAeda.sigma21HPD(:,5); s02]);
ylimL=min([outBAeda.sigma21HPD(:,4); s02]);
ylim([ylimL ylimU])
% Set xlim
xlim([xlimL xlimU]);
xL = get(my_subplot,'XLim');
% Add the horizontal line which corresponds to prior value of $\sigma^2$
line(xL,[s02 s02],'Color','r','LineWidth',lwd);
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
xlabel('Subset size m','FontSize',FontSize);
suplabel(['Housing data; forward plots in the interval ['...
num2str(xlimL) ',' num2str(xlimU) ...
']  of HPD regions for \beta and \sigma^2'],'t');
Observed curve of r_min is at least 10 times greater than 99.99% envelope
--------------------------------------------------
-------------------------
Signal detection loop
Tentative signal in central part of the search: step m=511 because
rmin(511,546)>99.99% and rmin(510,546)>99.99% and rmin(512,546)>99.99%
rmin(511,546)>99.999%

-------------------
Signal validation exceedance of upper envelopes
Validated signal
-------------------------------
Start resuperimposing envelopes from step m=510
Superimposition stopped because r_{min}(512,529)>99.9% envelope
Subsample of 528 units is homogeneous
----------------------------
Final output
Number of units declared as outliers=18
Summary of the exceedances
1          99         999        9999       99999
10          64          48          38          29

m=100
m=200
m=300
m=400
m=500


### Plot posterior estimates of beta and sigma2 in the last steps.

Plot posterior estimates of beta and sigma2 in the interval (subset size) [250 n+1].

% In this example for the house price data we monitor the forward plots
% of the parameters of the linear model adding 95% and 99% HPD
% regions. The first 8 panels refer to the elements of $\beta_1$
% and the bottom right-hand panel refer to the estimate of sigma2.
% The horizontal lines correspond to prior values
% The vertical lines refer to the step prior to the introduction of the first outlier
% init = point to start monitoring diagnostics along the FS
% run routine FSRB in order to find the outliers automatically
% setup parameters
n=size(hprice,1);
y=hprice(:,1);
X=hprice(:,2:5);
bayes=struct;
n0=5;
bayes.n0=n0;
% set \beta components
beta0=0*ones(5,1);
beta0(2,1)=10;
beta0(3,1)=5000;
beta0(4,1)=10000;
beta0(5,1)=10000;
bayes.beta0=beta0;
% \tau
s02=1/4.0e-8;
tau0=1/s02;
bayes.tau0=tau0;
% R prior settings
R=2.4*eye(5);
R(2,2)=6e-7;
R(3,3)=.15;
R(4,4)=.6;
R(5,5)=.6;
R=inv(R);
bayes.R=R;
% Automatic outlier detection procedure.
outBA=FSRB(y,X,'bayes',bayes', 'plots',0);
dout=n-length(outBA.ListOut);
% init = initial point to start monitoring
init=250;
xlimL=init; % lower value fo xlim
xlimU=n+1;  % upper value of xlim
out=FSRBeda(y,X,'bayes',bayes,'init',init, 'conflev', [0.95 0.99]);
% Set font size, line width and line style
figure;
lwd=2.5;
FontSize=14;
linst={'-','--',':','-.','--',':'};
for j=1:5
my_subplot=subplot(3,2,j);
hold('on')
% plot 95% and 99% HPD  trajectories
plot(out.beta1(:,1),out.beta1HPD(:,1:2,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','b')
plot(out.beta1(:,1),out.beta1HPD(:,3:4,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% plot posterior estimate of beta1_j
plot(out.beta1(:,1),out.beta1(:,j+1)','LineStyle',linst{1},'LineWidth',lwd,'Color','k')
% Add the horizontal line which corresponds to prior values
xL = get(my_subplot,'XLim');
line(xL,[beta0(j) beta0(j)],'Color','r','LineWidth',lwd);
% Set ylim
ylimU=max([out.beta1HPD(:,4,j); beta0(j)]);
ylimL=min([out.beta1HPD(:,3,j); beta0(j)]);
ylim([ylimL ylimU])
% Set xlim
xlim([xlimL xlimU]);
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
ylabel(['$\hat{\beta_' num2str(j-1) '}$'],'Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
if j>4
xlabel('Subset size m','FontSize',FontSize);
end
end
% Subplot associated with the monitoring of sigma^2
subplot(3,2,6);
%figure()
hold('on')
% 99%
plot(out.sigma21HPD(:,1),out.sigma21HPD(:,4:5),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% 95%
plot(out.sigma21HPD(:,1),out.sigma21HPD(:,2:3),'LineStyle',linst{2},'LineWidth',lwd,'Color','b')
% Plot 1/tau1
plot(out.S21(:,1),1./out.S21(:,3),'LineWidth',lwd,'Color','k')
ylabel('$\hat{\sigma}^2$','Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
% Set ylim
ylimU=max([out.sigma21HPD(:,5); s02]);
ylimL=min([out.sigma21HPD(:,4); s02]);
ylim([ylimL ylimU])
% Set xlim
xlim([xlimL xlimU]);
xL = get(my_subplot,'XLim');
% Add the horizontal line which corresponds to prior value of $\sigma^2$
line(xL,[s02 s02],'Color','r','LineWidth',lwd);
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
xlabel('Subset size m','FontSize',FontSize);
suplabel(['Housing data; forward plots in the interval ['...
num2str(xlimL) ',' num2str(xlimU) ...
']  of HPD regions for \beta and \sigma^2'],'t');
Observed curve of r_min is at least 10 times greater than 99.99% envelope
--------------------------------------------------
-------------------------
Signal detection loop
Tentative signal in central part of the search: step m=511 because
rmin(511,546)>99.99% and rmin(510,546)>99.99% and rmin(512,546)>99.99%
rmin(511,546)>99.999%

-------------------
Signal validation exceedance of upper envelopes
Validated signal
-------------------------------
Start resuperimposing envelopes from step m=510
Superimposition stopped because r_{min}(512,529)>99.9% envelope
Subsample of 528 units is homogeneous
----------------------------
Final output
Number of units declared as outliers=18
Summary of the exceedances
1          99         999        9999       99999
10          64          48          38          29

m=100
m=200
m=300
m=400
m=500


### Plot of HPDI for BankProfit data.

XX=load('BankProfit.txt');
X=XX(:,1:end-1);
y=XX(:,end);
beta0=zeros(10,1);
beta0(1,1)=-0.5;        %
beta0(2,1)=9.1;         % Number of products (NUMPRO)
beta0(3,1)=0.001;       % direct revenues  (DIRREV)
beta0(4,1)=0.0002;      % indirect revenues   (INDREV)
beta0(5,1)=0.002;       %  savings accounts   SAVACC
beta0(6,1)=0.12;        %  number of operations   NUMOPE
beta0(7,1)=0.0004;      %  total amount of operations  TOTOPE
beta0(8,1)=-0.0004;     %  Bancomat POS
beta0(9,1)=1.3;         %  Number of cards   NUMCAR
beta0(10,1)=0.00004;    %  Amount in cards   TOTCAR
% \tau
s02=10000;
tau0=1/s02;
% number of obs in which prior was based
n0=1500;
bayes=struct;
bayes.R=R;
bayes.n0=n0;
bayes.beta0=beta0;
bayes.tau0=tau0;
n=length(y);
% run routine FSRB in order to find the outliers automatically
outBA=FSRB(y,X,'bayes',bayes', 'plots',0);
dout=n-length(outBA.ListOut);
% init = point to start monitoring diagnostics along the FS
init=1700;
xlimL=init;
xlimU=n+1;
outBAeda=FSRBeda(y,X,'bayes',bayes,'init',init, 'conflev', [0.95 0.99]);
% Set font size, line width and line style
figure;
lwd=2.5;
FontSize=14;
linst={'-','--',':','-.','--',':'};
nr=4;
nc=3;
for j=1:length(beta0)
my_subplot=subplot(nr,nc,j);
hold('on')
% plot 95% and 99% HPD  trajectories
plot(outBAeda.beta1(:,1),outBAeda.beta1HPD(:,1:2,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','b')
plot(outBAeda.beta1(:,1),outBAeda.beta1HPD(:,3:4,j),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% plot posterior estimate of beta1_j
plot(outBAeda.beta1(:,1),outBAeda.beta1(:,j+1)','LineStyle',linst{1},'LineWidth',lwd,'Color','k')
% Add the horizontal line which corresponds to prior values
xL = get(my_subplot,'XLim');
line(xL,[beta0(j) beta0(j)],'Color','r','LineWidth',lwd);
% Set ylim
ylimU=max([outBAeda.beta1HPD(:,4,j); beta0(j)]);
ylimL=min([outBAeda.beta1HPD(:,3,j); beta0(j)]);
ylim([ylimL ylimU])
% Set xlim
xlim([xlimL xlimU]);
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
ylabel(['$\hat{\beta_' num2str(j-1) '}$'],'Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
if j>nr*(nc-1)
xlabel('Subset size m','FontSize',FontSize);
end
end
% Subplot associated with the monitoring of sigma^2
subplot(4,3,11);
%figure()
hold('on')
% 99%
plot(outBAeda.sigma21HPD(:,1),outBAeda.sigma21HPD(:,4:5),'LineStyle',linst{4},'LineWidth',lwd,'Color','r')
% 95%
plot(outBAeda.sigma21HPD(:,1),outBAeda.sigma21HPD(:,2:3),'LineStyle',linst{2},'LineWidth',lwd,'Color','b')
% Plot 1\/tau1
plot(outBAeda.S21(:,1),1./outBAeda.S21(:,3),'LineWidth',lwd,'Color','k')
ylabel('$\hat{\sigma}^2$','Interpreter','LaTeX','FontSize',20,'rot',-360);
set(gca,'FontSize',FontSize);
% Set ylim
ylimU=max([outBAeda.sigma21HPD(:,5); s02]);
ylimL=min([outBAeda.sigma21HPD(:,4); s02]);
ylim([ylimL ylimU])
ylimL=8000;
ylimU=14000;
ylim([ylimL ylimU])
% Set xlim
xlim([xlimL xlimU]);
xL = get(my_subplot,'XLim');
% Add the horizontal line which corresponds to prior value of $\sigma^2$
line(xL,[s02 s02],'Color','r','LineWidth',lwd);
% Add vertical line in correspondence of the step prior to the
% entry of the first outlier
line([dout; dout],[ylimL; ylimU],'Color','r','LineWidth',lwd);
xlabel('Subset size m','FontSize',FontSize);
suplabel(['Bank profit data; forward plots in the interval ['...
num2str(xlimL) ',' num2str(xlimU) ...
']  of HPD regions for \beta and \sigma^2'],'t');

## Input Arguments

### y — A vector with n elements that contains the response variables. Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) with missing or infinite values will automatically be excluded from the computations.

Data Types: single| double

### X — Data matrix of explanatory variables (also called 'regressors') of dimension (n x p-1). Rows of X represent observations, and columns represent variables.

Missing values (NaN's) and infinite values (Inf's) are allowed, since observations (rows) with missing or infinite values will automatically be excluded from the computations.

PRIOR INFORMATION $\beta$ is assumed to have a normal distribution with mean $\beta_0$ and (conditional on $\tau_0$) covariance $(1/\tau_0) (X_0'X_0)^{-1}$.

$\beta \sim N( \beta_0, (1/\tau_0) (X_0'X_0)^{-1} )$

Data Types: single| double

### Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as  Name1,Value1,...,NameN,ValueN.

Example:  'intercept',false , bayes=struct;bayes.R=R;bayes.n0=n0;bayes.beta0=beta0;bayes.tau0=tau0; , bsb=[2 5 1]; , 'init',100 starts monitoring from step m=100 , 'nocheck',true , 'conflev',[0.90 0.93] 

### intercept —Indicator for constant term.true (default) | false.

Indicator for the constant term (intercept) in the fit, specified as the comma-separated pair consisting of 'Intercept' and either true to include or false to remove the constant term from the model.

Example:  'intercept',false 

Data Types: boolean

### bayes —It specifies prior information.structure.

A structure which specifies prior information.

If structure bayes is not supplied the default values which are used are: beta0= zeros(p,1) (vector of zeros); R=eye(p) (Identity matrix); tau0=1/1e+6 (very large value for the prior variance, that is a very small value for tau0); n0=1 (just one prior observation).

Strucure bayes contains the following fields:

Value Description
beta0

p-times-1 vector containing prior mean of $\beta$

R

p-times-p positive definite matrix which can be interepreted as $X_0'X_0$ where $X_0$ is a n0 x p matrix coming from previous experiments (assuming that the intercept is included in the model)

tau0

scalar. Prior estimate of $\tau=1/ \sigma^2=a_0/b_0$.

The prior distribution of tau0 is a gamma distribution with parameters a and b, that is $p(\tau_0) \propto \tau^{a_0-1} \exp (-b_0 \tau)$;

$E(\tau_0)= a_0/b_0$.

n0

scalar. Sometimes it helps to think of the prior information as coming from n0 previous experiments.

Therefore we assume that matrix $X_0$ (which defines R), was made up of n0 observations.

Example:  bayes=struct;bayes.R=R;bayes.n0=n0;bayes.beta0=beta0;bayes.tau0=tau0; 

Data Types: double

### bsb —list of units forming the initial subset.vector.

if bsb=0 then the procedure starts with p units randomly chosen else if bsb is not 0 the search will start with m0=length(bsb). The default value of bsb is '' that is in the first step just prior information is used.

Example:  bsb=[2 5 1]; 

Data Types: double

### init —Search initialization.scalar.

scalar, specifies the point where to start monitoring required diagnostics. if init is not specified it will be set equal to :

p+1, if the sample size is smaller than 40;

min(3*p+1,floor(0.5*(n+p+1))), otherwise.

Example:  'init',100 starts monitoring from step m=100 

Data Types: double

### nocheck —Check input arguments.boolean.

Scalar. If nocheck is equal to true no check is performed on matrix y and matrix X. Notice that y and X are left unchanged. In other words the additional column of ones for the intercept is not added. As default nocheck=false. See routine chkinputRB for the details of the operations.

Example:  'nocheck',true 

Data Types: boolean

### conflev —confidence levels to be used to compute HPDI.vector.

This input option is used just if input stats=1. The default value of conflev is [0.95 0.99] that is 95% and 99% HPDI confidence intervals are computed.

Example:  'conflev',[0.90 0.93] 

Data Types: double

## Output Arguments

### out — description Structure

T consists of a structure 'out' containing the following fields:

Value Description
RES

Scaled residuals. Matrix.

n x (n-init+1) matrix containing the monitoring of scaled residuals.

1st row = residual for first unit;

...;

nth row = residual for nth unit.

LEV

Leverage. Matrix.

(n+1) x (n-init+1) = matrix containing the monitoring of leverage.

1st row = leverage for first unit;

...;

nth row = leverage for nth unit.

BB

n x (n-init+1) matrix containing the information about the units belonging to the subset at each step of the forward search.

1st col = indexes of the units forming subset in the initial step;

...;

last column = units forming subset in the final step (all units).

mdrB

n-init x 3 matrix which contains the monitoring of Bayesian minimum deletion residual or (m+1)ordered residual at each step of the forward search.

1st col = fwd search index (from init to n-1);

2nd col = minimum deletion residual;

3rd col = (m+1)-ordered residual.

Remark: these quantities are stored with sign, that is the min deletion residual is stored with negative sign if it corresponds to a negative residual.

msrB

n-init+1 x 32= matrix which contains the monitoring of maximum studentized residual among units belonging to subset.

1st col = fwd search index (from init to n);

2nd col = maximum studentized residual.

Coo

(n-init+1) x 2 matrix containing the monitoring of Cook or modified Cook distance in each step of the forward search.

1st col = fwd search index (from init to n);

2nd col = monitoring of Cook distance.

beta1

(n-init+1) x (p+1) matrix containing the monitoring of posterior mean of $\beta$ (regression coefficents) $\beta_1 = (c*R + X'X)^{-1} (c*R*\beta_0 + X'y)$

Gam

(n-init+1) x 3 matrix containing:

1st col = fwd search index (from init to n);

2nd col = shape parameter $a_1$ of the posterior gamma distribution of tau;

3rd col = scale parameter $b_1$ of the posterior gamma distribution of tau.

Remark: $a_1 = 0.5 (c*n0 + m)$ where m is subset size;

$b_1 = 0.5 * ( n0 / \tau_0 + (y-X*\beta_1)'y +(\beta_0-\beta_1)'*c*R*\beta_0 )$

covbeta1

p x p x (n-init+1) 3D array containing posterior covariance matrix (conditional on $tau_1$) of $\beta$.

$cov(\beta_1) = (1/tau_1) * (c*R + X'X)^{-1}$;

where $tau_1$ is defined as $a_1/b_1$ (that is through the gamma parameters of the posterior distribution of $\tau$).

The posterior distribution of $\tau$ is a gamma distribution with parameters $a_1$ and $b_1$.

S21

(n-init+1) x 3 matrix containing the monitoring of posterior estimate of $\sigma^2$ and $\tau$ in each step of the forward search.

1st col = fwd search index (from init to n);

2nd col = monitoring of $\sigma^2_1$ (posterior estimate of $\sigma^2$);

3rd col = monitoring $\tau_1$ (posterior estimate of $\tau$).

Bpval

(n-init+1) x (p+1) containing Bayesian p-values.

p-value = $P(|t| > | \hat{\beta} / se(beta) |)$ = prob. of beta different from 0.

1st col = fwd search index (from init to n);

2nd col = p-value for first variable;

...;

(p+1) col = p-value for p-th variable.

beta1HPD

(n-init+1)-by-2-by-p 3D array.

Bhpd(:,:,1) = lower and upper HPDI of first element of $\beta_1$ (posterior estimate of $\beta$);

...;

Bhpd(:,:,p) = lower and upper HPDI of last element of $\beta_1$ (posterior estimate of $\beta$).

tau1HPD

(n-init+1) x 3 containing HPDI for $\tau_1$.

1st col = fwd search index (from init to n);

2nd col = lower value of HPDI;

3rd col = upper value of HPDI;

sigma21HPD

(n-init+1) x 3 containing HPDI for $\sigma^2_1$.

1st col = fwd search index (from init to n);

2nd col = lower value of HPDI;

3rd col = upper value of HPDI.

postodds

(n-init+1)-by-(p+1) matrix which contains posterior odds for $\beta_j=0$.

For example the posterior odd of $\beta_0=0$ is p(y| model which contains all expl variables except the one associated with beta0) divided by p(y| model which contains all expl variables).

1st col = fwd search index (from init to n);

2nd col = posterior odd for $beta_1$;

...;

(p+1) col = posterior odd for $beta_p$.

modelprob

(n-init+1)-by-(p+1) matrix which contains posterior model probability of the model which excludes variable j. For example if modelprob(j)= 0.28, that is if the probability of the model which does not contain variable j is equal to 0.28, it means that there is a 28 chance that beta_j=0 and a 72 chance that it is not.

1st col = fwd search index (from init to n);

2nd col = posterior model prob of the model which excludes $\beta_1$;

...;

(p+1) col = posterior model prob of the model which excludes $beta_p$.

Un

(n-init) x 11 matrix which contains the unit(s) included in the subset at each step of the fwd search.

REMARK: in every step the new subset is compared with the old subset. Un contains the unit(s) present in the new subset but not in the old one Un(1,2) for example contains the unit included in step init+1 Un(end,2) contains the units included in the final step of the search.

y

A vector with n elements that contains the response variable which has been used.

X

Data matrix of explanatory variables which has been used (it also contains the column of ones if input option intercept was missing or equal to true).

class

'FSRBeda'.

## References

Chaloner, K. and Brant, R. (1988), A Bayesian Approach to Outlier Detection and Residual Analysis, "Biometrika", Vol. 75, pp. 651-659.

Riani, M., Corbellini, A. and Atkinson, A.C. (2018), Very Robust Bayesian Regression for Fraud Detection, "International Statistical Review", http://dx.doi.org/10.1111/insr.12247

Atkinson, A.C., Corbellini, A. and Riani, M., (2017), Robust Bayesian Regression with the Forward Search: Theory and Data Analysis, "Test", Vol. 26, pp. 869-886, DOI 10.1007/s11749-017-0542-6