Author

Keywords

Review Status

1 Introduction

This tutorial explains how to implement a Particle Swarm (PS) optimisation for robust regression as used in Gilli and Schumann [2]. It is still ongoing research, so comments are highly welcome.

The chosen examples are LTS- and LMS-regression.

2 The algorithm

The PS algorithm is implemented as a function PSO that is called with three arguments:

PSO(pso,dataList,OF)

pso is a list that holds the settings of the PS algorithm. dataList is a list that holds the pieces of data necessary to evaluate the objective function. For a regression model, dataList will typically contain the matrix $X$ and the vector $y$. OF is the objective function. In principle, we need not pass these data, since the scoping rules in R also allow to access objects outside a function's environment. However, it makes for a clearer structure of the algorithm to pass all required objects explicitly to the function.

The function will return a list of two elements: beta which is the estimated parameter vector of the model, and OFvalue, the objective function value associated with this solution.

Settings

PS requires us to set several parameters. These are collected in the list pso and then passed to the function. For a discussion of these parameters, see Gilli and Schumann [2]. Initial velocity initV serves as a factor to the randomly set the velocities (computed as $\mathtt{initV} \times \mathtt{mRU(d,nP)}$).

The objective function

Our solutions are stored in a matrix mP which we pass to the objective function. Alternatively, we could have called a function of the apply-family from the main algorithm. However, for objective functions that could be vectorised, this setup here allows for easier reuse of the code. Here, we use the Least Trimmed Squares (LTS) estimator.

In order not to depend on a single data set, we ran 50 times the data-generating model of Salibian-Barrera and Yohai [4] (see the function genData given above) with 400 observations and 40 outliers (set to 150). For each data set, we estimated coefficients with ltsReg and PSO and compute the quantity

for the residuals from ps. The $h$th order statistic of the squared residuals is defined by

h <- (nO+nP+1) %/% 2

The following figure gives the distribution of the ratio computed in Equation X for a model with three regressors. A value of 1 implies both models gave the same answer in terms of their objective functions, a value below unity indicates that the PS algorithm performed better.

Another criterion

Let us demonstrate the flexibility of PS by replacing the objective function: we now use the Least Median of Squares (LMS) estimator. The only thing we need to do is to rewrite the objective function.