PF_lm_ss: Parameter Estimation Of Linear Regression Using Particle Filters With Simple Sampling

Description

Estimation of the coefficients of a linear regression based on the particle filters algorithm. This function is similar to PF_lm except for the resampling method which in this case is the simple sampling. As a result, the user can try higher number of particles.

Usage

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

Value

A list containing the following elemnts:

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.

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 variance of the normal likelihood, 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

Author

Christian Llano Robayo, Nazrul Shaikh.

Details

Estimation of the coefficients of a linear regression:

$Y = β_{0} + β_{1} * X_{1} + . . . + β_{n} * X_{n} + ϵ$ , ( $ϵ \sim N (0, 1)$ )

using particle filter methods. The state-space equations are:

$(E q . 1) X_{k} = a_{0} + a_{1} * X 1_{k - 1} + . . . + a_{n} * X n_{k - 1}$

$(E q . 2) Y_{k} = X_{k} + ϵ, ϵ \sim N (0, σ)$

where, k = 2, ... , number of observations; $a_{0}, . . ., a_{n}$ are the parameters to be estimated (coefficients), and $σ = 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 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 $ϵ$ together with the coefficients. If initDisPar is missing, the initial priors are taken using lm() and coeff() plus-minus lbd as a reference.

The resampling method used corresponds to the simple sampling, i.e., we take a sample of n particles with probability equals the likelihood computed on each iteration. In addition, on each iteration white noise is added to avoid particles to degenerate.

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} + ϵ$

References

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

Run this code


if (FALSE) {
#### Using default Data1, no sigma estimation ####
Res <- PF_lm_ss(n = 10000L, sigma_est = FALSE) #10 times more than in PF_lm
lapply(Res,class) # Structure of returning list.
###Summary of estimated parameters
Res$summ
#Evolution of the estimated parameters
par(mfrow=c(2, 2))
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_ss(n = 1000L, sigma_est = TRUE)
lapply(Res2,class) # Structure of returning list
###Summary of the estimated parameters
Res2$summ
#Evolution of the estimated parameters
par(mfrow=c(2, 3))
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_ss(n = 10000L, sigma_est = FALSE, initDisPar = b0)
lapply(Res3,class) # Structure of returning list.
###Summary of the estimated parameters
Res3$summ
#Evolution of the estimated parameters
par(mfrow=c(2, 2))
for (i in 1:4){
plot(apply(Res3$stateP_res[[i]],1,mean), main = colnames(Res3$summ)[i], col="blue",
    xlab =  "", ylab = "",type="l")
    }
}

Run the code above in your browser using DataLab

Get 50% off unlimited learning