PF_lm {LMfilteR}

Parameter Estimation Of Linear Regression Using Particle Filters

Description

Estimation of the coefficients of a linear regression based on the particle filters algorithm (PF); given Data1, the function estimates the coefficients of the regression as the samples (particles) of coefficients that are more likely to occur. Optional, the standard deviation of the error term can also be estimated. Synthetic data is generated in case no provided Data1.

Usage

PF_lm(Y, Data1, n = 500L, sigma_l = 1, sigma_est = FALSE, initDisPar, lbd = 2)

Arguments

Y

numeric. The response variable

Data1

matrix. The matrix containing the independent variables

n

integer. Number of particles, by default 500

sigma_l

The value of the variance in the likelihood normal, 1 by default

sigma_est

logical. If TRUE takes the last row of initDisPar as prior estimation of the standard deviation, see more in Details

initDisPar

matrix. Values a, b of the uniform distribution (via runif) for each parameter to be estimated, see more in Details

lbd

numeric. A number to be added and substracted from the priors when initDisPar is not provided

Details

Estimation of the coefficients of a linear regression:

Y = \beta_0 + \beta_1*X_1 + ... + \beta_n*X_n + \epsilon, (\epsilon \sim N(0,1))

using particle filter methods. The implementation follows An, D., Choi, J. H., Kim, N. H. (2013) adapted for the case of linear models.

The state-space equations are:

(Eq. 1) X_{k} = a_0 + a_1 * X1_{k-1} + ... + a_n * Xn_{k-1}

(Eq. 2) Y_{k} = X_{k} + \epsilon, \epsilon \sim N(0,\sigma)

where, k = 2, ... , number of observations; a_0, ... , a_n are the parameters to be estimated (coefficients), and \sigma = 1 .

The priors of the parameters are assumed uniformly distributed. initDisPar is a matrix for which the number of rows is the number of independent variables plus one when sigma_est = FALSE, (plus one since we also estimate the constant term of the regression), or plus two when sigma_est = TRUE (one for the constant term of the regression, and one for the estimation of sigma). The first and second column of initDisPar are the corresponding arguments a and b of the uniform distribution (stats::runif) of each parameter prior. The first row of initDisPar is the prior guess of the constant term. The following rows are the prior guesses of the coefficients. When sigma is estimated, i.e., if sigma_est = TRUE the last row of initDisPar corresponds to the prior guess for the standard deviation. If sigma_est = FALSE, then the standard deviation of the likelihood is sigma_l. If sigma_est = TRUE, the algorithm estimates the standard deviation of \epsilon together with the coefficients. If initDisPar is missing, the initial priors are taken using lm() and coeff() plus-minus lbd as a reference.

The cumulative density function method is used for resampling.

In case no Data1 is provided, a synthetic data set is generated automatically taking three normal i.i.d. variables, and the dependent variable is computed as in Y = 2 + 1.25*X_1 + 2.6*X_2 - 0.7*X_3 + \epsilon

Value

A list with the following elements:

stateP_res: A list of matrices with the PF estimation of the parameters on each observation; the number of rows is the number of observations in Data1 and the number of columns is n.

Likel: A matrix with the likelihood of each particle obtained on each observation.

CDF: A matrix with the cumulative distribution function for each column of Likel.

summ: A summary matrix with statistics of the estimated parameters.

Author(s)

Christian Llano Robayo, Nazrul Shaikh.

References

An, D., Choi, J. H., Kim, N. H. (2013). Prognostics 101: A tutorial for particle filter-based prognostics algorithm using Matlab. Reliability Engineering & System Safety, 115, DOI: https://doi.org/10.1016/j.ress.2013.02.019

Ristic, B., Arulampalam, S., Gordon, N. (2004). Beyond the Kalman filter: particle filters for tracking applications. Boston, MA: Artech House. ISBN: 158053631X.

West, M., Harrison, J. (1997). Bayesian forecasting and dynamic models (2nd ed.). New York: Springer. ISBN: 0387947256.

Examples

## Not run: 
#### Using default Data1, no sigma estimation ####
Res <- PF_lm(n=1000L, sigma_est = FALSE)
lapply(Res,class) # Structure of returning list.
###Summary of estimated parameters
Res$summ
par(mfrow=c(2, 2)) #Evolution of the estimated parameters
for (i in 1:4){
plot(apply(Res$stateP_res[[i]],1,mean), main = colnames(Res$summ)[i], col="blue",
      xlab =  "", ylab = "",type="l")
}


#### Using default Data1, with sigma estimation ####
Res2 <- PF_lm(n=1000L, sigma_est = TRUE)
lapply(Res2,class) # Structure of returning list
sumRes2 <- sapply(2:6, function(i)
###Summary of estimated parameters
Res2$summ

par(mfrow=c(2, 3)) #Evolution of the estimated parameters
for (i in 1:5){
plot(apply(Res2$stateP_res[[i]],1,mean), main = colnames(Res2$summ)[i], col="blue",
       xlab =  "", ylab = "",type="l")
}


#### Using default Data1, given initDisPar ####
b0 <- matrix(c( 1.9, 2, # Prior of a_0
                1, 1.5, # Prior of a_1
                2, 3,   # Prior of a_2
                -1, 0)  # Prior of a_3
                ,ncol = 2, byrow = TRUE )
Res3 <- PF_lm(n=500L, sigma_est = FALSE, initDisPar = b0)
lapply(Res3,class) # Structure of returning list.

###Summary of estimated parameters
Res3$summ
par(mfrow=c(2, 2)) #Evolution of the estimated parameters
for (i in 1:4){
plot(apply(Res3$stateP_res[[i]],1,mean), main = colnames(Res3$summ)[i], col="blue",
      xlab =  "", ylab = "",type="l")
}

## End(Not run)