1 Introduction
In its classical form, the multiarmed bandit (MAB) problem requires that a learning agent makes a choice, or selects an action, from alternatives, across
trials, for each trial. After a choice is made in a given trial, the system presents the learning agent with a numerical reward from a stationary probability distribution associated with the action taken. The goal of the learning agent is to dynamically and sequentially choose alternatives, referred to as arms, which will maximize the expected total reward across the
trials [Sutton and Barto1998]. An extended version of this problem statement is the linear bandit problem, which is sometimes referred through other names such as linear stochastic bandits or linear parameterized bandits, and sometimes studied under specific settings such as contextual or associative bandits with linear payoff functions. In its general form, originally proposed by Auer 2 (whereas variants of the problem was considered prior to that [Woodroofe1979, Kaelbling1994]), each arm or alternative is first parameterized into a set of features, which is known to the learning agent. The reward is still a sample from a stationary probability distribution, but its expected value is seen as the inner product of the feature vector and a fixed weight vector, which indicates the influence of each feature on the expected reward. The learning agent therefore seeks to understand these weights in its exploration phase, and uses this to exploit the system by picking promising arms. The critical advantage of this approach is that since we are tagging the rewards to the features, and not the arms directly, we do not need to pull each arm (the set of which could be large or even infinite), but would still be able to learn about the expected reward of each arm by understanding them through a common set of features (which is expected be smaller in size).
[Auer2002, Dani et al.2008, Rusmevichientong and Tsitsiklis2010, AbbasiYadkori et al.2011, Li et al.2010].The linear stochastic bandit conception allows us to extend various realworld online learning problems to a bandit framework. Studies by Abe and Long 7 , Auer 2, and Li et al., 6 discuss the application of a selection problem in displaying internet banner Ads or News articles. The agent needs to choose an Ad among many, to display to a specific user. Here the agent focuses its learning on understanding the features, which are combination of the User’s and Ad’s characteristics, as opposed to learning directly about each individual ad. The reward is received when a user clicks on the Ad. Similarly, Rusmevichientong and Tsitsikilis 4 presents an application in marketing where the agent is tasked with choosing a product (arm) to offer to a customer, and the various product characteristics, such as price and popularity, are the features. The scope for application also goes beyond selection problems that span single objects or single processes conceptualized as arms (products, Ads). Situations exist where we have a set of decisions, or a combination of actions that need to be conceptually repackaged as a single arm. As discussed in Dani et al., 3, take the typical case of choosing one out of clinical treatments, which is modelled as a classical bandit problem with arms. If we wanted to extend this context to a decision problem where we could choose any combination of the clinical treatments to be applied on a subject, then the set of possible decisions increase to (as there can be interactions between the treatments). In this situation, we could treat each of the combinations as separate arms in a linear bandit framework. We could conceptualize, at minimum, features (corresponding to presence of absence of a treatment), and perhaps even some interaction terms and higher order transformations as additional features. Another novel application related to such combinatorics is discussed by Awerbuch and Kleinberg 8 where an agent needs to figure out, online, the quickest path of getting from one point to another in a network. Here, multiple paths involving a combination of edges in the network can be selected to achieve this task and the agent learns the cost of choosing a certain path, but not the cost of each of the edges. Modelling this through linear bandit framework involves considering each path as an arm, and each edge as a feature.
In this study, we propose two algorithms to solve the linear stochastic bandit problems without making any assumptions about the noise in the system. Inspired by practices in the statistics community to generate confidence bands on linear models through bootstraps [Thompson1978], we apply these bootstrap algorithms in conjunction with our linear bandit formulation to create an arm pulling agent. Both these algorithms rely on the popular concept of using an upper confidence bound (UCB) associated with arms as described in Auer 2. The UCBs are derived from the bootstraps and used to choose an arm. The algorithms can also be conceptualized as an alternate way of building confidence sets as described in AbbasiYadkori et al. 5. We compare these algorithms to other well researched solutions in the linear bandit space. Specifically, we look at the algorithm OFUL [AbbasiYadkori et al.2011], LinUCB [Li et al.2010], and Thompson sampling for linear bandits [Agrawal and Goyal2012]. We perform this comparison using a simulation on a hierarchical probability metamodel, which create response surfaces inspired by real world systems.
As described later in the Section 3, this approach allows us to choose the nature of the noise in the system. We experimented with different noise models, paying special attention to a case where the distributional assumptions of the stateoftheart methods are violated. When we use Gaussian noise, which fit the assumptions made by the methods, the performance of our approach is comparable to the stateoftheart. When we move away from this to something that is not subGaussian (we chose the Laplace distribution in our study) the proposed methods start outperforming the stateoftheart.
The contributions of this study are as follows:
1) This is the first ever use of bootstrap estimates of the confidence bounds in a bandit setup. Earlier contributions are limited to estimation of a posterior distribution or evaluation of solutions
[Eckles and Kaptein2014, Osband and Roy2015]. The chief advantage of using bootstraps is that we make no distributional assumption about the stochasticity in the environment. We demonstrate the advantage of this by showing empirically that our approaches out perform stateoftheart methods when the underlying assumptions are violated.2) While problems with combinatorial sets of decisions have been modelled as linear bandits earlier, the systems have been limited to modelling additive effects. We present a systematic method for including higher order interactions between choices. This allows us to expand the applicability of bandits to domains such as design of experiments where models with such higher order effects are necessary.
3) The other contribution in this study pertains to the test environment used for the algorithms. Similar to other empirical studies on bandits, this study also uses synthetic data. However, the simulation data that is generated by our system is derived from realworld data gathered from published experiments on engineering systems. This metamodelling approach allows us to capture idiosyncrasies of the realworld environments, but also perform simulations on a large number of response surfaces (as opposed to case studies). This study seeks to demonstrate the use of such an approach to testing bandit algorithms and also provide users with a specific response surface generating environment.
The paper is organized as follows. Section 2 introduces the proposed bootstrap algorithms, the notation and pseudo code, and the comparison algorithms. Section 3 describes the simulation environment. Section 4 presents the results with a discussion and identifies future research directions.
2 Algorithms
2.1 Discussion
Bootstrapping as a statistical approach seeks to sample from the data with replacement. This creates a version of the data which is different from the original set and therefore exploits the central idea of bootstrapping that The population is to the sample, as the sample is to the bootstrap sample’. The two main algorithms presented in this study are the Fixed and
random Bootstrap Bandits, corresponding to the similarly named approaches of bootstrapping to evaluate the confidence bands (different from confidence bounds as it applies to a model not a single parameter) of linear regression models discussed in the statistics literature. The use of these bootstrap approaches, especially with the above mentioned terminology, is first seen in Thompson 10 and well discussed in various other sources
[Efron and Tibshirani1994, Breiman1992, Breiman and Spector1992].The algorithms presented in this study build on Auer 2’s idea of using the upper confidence bound (UCB) to create an arm pulling agent. The upper confidence bound is constructed using the bootstrap approach on linear regression as discussed above. The proposed algorithm can also be seen as an alternate way of implementing OFUL [AbbasiYadkori et al.2011], which in turn also builds on the UCB idea of Auer2. In trial the OFUL algorithm seeks to build the confidence set for the parameters or weight vector , based off of rewards and arms pulled in previous rounds up to . This set then serves as a constraint in the maximization of the inner product, where and , where is the set of all unique arms. The arm resulting from this maximization is then selected or pulled for round . The confidence set is built by constructing an ellipsoid in the parameter space around the regularized leastsquares estimate of , such that with a high predefined probability lies in . Our bootstrap approach does not create this ellipsoid but seeks to directly populate or create a confidence set by using the bootstraps to create multiple linear regression fits, and therefore multiple sets of parameters. This study adopts the same optimisminthefaceofuncertainty optimisation to choose the arm .
In a bandit framework, the use of sampling directly from past data to create an arm pulling agent is rarely seen (Thompson sampling, for instance, samples from a distribution, which is created from the data). However, there are some notable exceptions [Baransi et al.2014, Eckles and Kaptein2014, Osband and Roy2015]. In Baranasi et al., 17 a form of subsampling without replacement is applied to the MAB problem. Whereas, the technical notes by Eckles and Kaptein 24, and Osband and Roy 25 look at computationally efficient modifications to the core bootstrap idea and apply it to generate posterior distributions in Thompson sampling framework for the MAB problem. There are various advantages to using the bootstrap approaches. The main advantage being that these methods allow us to get confidence bounds without making any assumptions on the distributions or independence of variates, unlike some of the baseline methods. While such an approach could be computationally intensive, it could be very helpful in cases where the sampling distribution is unknown, or difficult to derive.
Similar to most bandit algorithms, the proposed solutions also require an initialization step. However, owing to the linear bandit conception, the algorithms described in this section can start by pulling a relatively small subset of the set of all unique arms. The number and selection of arms is based on the number of features the linear bandit algorithm seeks to estimate. At minimum the algorithm would need to pull more arms than the number of features. Also, the arms need to be carefully selected in order to ensure that the correlation between features (mulitcollinearity) is zero. This problem is well handled for discrete or discretized input spaces in offline design of experiments. A simple class of designs (sets of arms) that can satisfy the constraints discussed are Orthogonal Arrays. These are matrices which, in addition to satisfying more restrictive mathematical properties, select sets of arms were each feature assumes each discrete state an equal number of times and the correlation between any two features is zero [Montogomery1984]. In this study we propose the initialization step to be handled through orthogonal arrays of minimum size such that all parameters can be estimated.
2.2 Notation and Steps
Assume a system where there are unique arms, and features that describe each of these arms. In this context we are interested in two different sets of arms, represented by two different sets of vectors. We define the first vector set as which is of size and the comprehensive set of all unique arms. Each vector is of length corresponding to the number of features. We can also define the vector set , in matrix form, which represents the arms that have been tried, and for which rewards have been gathered. Matrix is of dimensions . Each row corresponds to the arm pulled at a particular trial or round, and after trials this matrix has rows and each arm is represented through its features. We present the vector of rewards as vector which is of also of size and corresponds to the rewards received from the each row of matrix . In addition to this we use the parameter which is tuneable and describes the percentage level of the upper confidence bound, and therefore implicitly the degree of exploration versus exploitation.
In both algorithms the step of initialization by pulling the minimum, carefully selected arms to estimate the features is conducted in step 2. In algorithm 1, the key steps for bootstrapping are captured in step 5 through 7. In algorithm 2, the key steps for bootstrapping are captured in steps 7 and 8. The next two paragraphs discuss the algorithms in detail.
In Algorithm 1, Random, for each trial we create a bootstrap sample of the data. Here, each input(arm)output(reward) pair is selected randomly with replacement, creating a bootstrapped dataset of the same original size (steps 56). This dataset is regressed to estimate the parameters corresponding to that bootstrap (step 7). This creates a bootstrapped sample of parameters. For this sample we estimate the expected reward across all arms (loop in steps 810). We repeat these steps for various bootstraps of the residuals (loop in steps 411).The principle of upper confidence is then applied by tagging each arm with its ”optimistic” performance (step 12). The arm with the best resulting reward is selected in that trial and receives a reward (step 14). The updated input matrix and reward vector are used for subsequent trials (step 3). The algorithm ends after exhausting all trials (step 16).
In the case of Algorithm 2, fixed, for each trial we fit a linear model to the historic data (step 4). The algorithm then computes a vector of residuals resulting from the difference between actual reward and the rewards expected from the linear model (step 5). The residuals are then bootstrapped to create new sample of residuals which are of the same size as the original set (step 7).Note that this is different from algorithm 1 which bootstraps the data points directly. We then apply the bootstrapped residuals to the fitted model to create new outputs, and reregress the new outputs to the fixed inputs (step 8). The algorithm then proceeds identically to the Random algorithm.
Literature suggests that the conceptualization of ‘what constitutes the sample’ should decide which of the two approaches we could use. If we believe that the inputs are fixed, and the output is a random sample from a distribution whose mean is determined by , then the Fixed is appropriate. However, if we believe that the inputs characterized by the matrix is itself a random sample, then we would find the Random to be appropriate [Efron and Tibshirani1994]. Given that in our context the initialization step involves a carefully selected input matrix in the form of an orthogonal array, Fixed should be a better choice, initially. However, subsequent arm pulls are based off of the algorithms reacting to reward data which is conceptualized as a random sample taken from s stationary probability distribution and hence Random would seem to be an appropriate choice, asymptotically.
2.3 Pseudo code
2.4 Baseline Algorithms
This study compares the bootstrapped linear bandits to three other established methods of working with linear bandits. These include OFUL [AbbasiYadkori et al.2011], LinUCB [Li et al.2010], and Thompson sampling for linear bandits [Agrawal and Goyal2012]. The OFUL algorithm, owing to its extensibility to the proposed algorithms, has been discussed in detail in section 2.1. The LinUCB algorithm uses a similar conception as OFUL and other UCB based algorithms of using the uncertainty in the regularized least squares linear model that is built from past data to estimate the upper confidence bound. The Thompson sampling approach utilizes a Bayesian set up, where, in each step a parameter vector is sampled from posterior distributions of the parameters. The arm that maximizes the reward for this sampled vector is chosen, and a corresponding reward is received. The posterior distribution is then updated to account for the newly received reward.
3 Testing Environment
Inspired by the examples discussed in [Dani et al.2008] and [Awerbuch and Kleinberg2004], we create a combinatorial application of linear bandits, and perform a simulation study on it. Specifically we replicate the extension of experimenting with clinical treatments. We take a sample case where treatments can be offered to patients, but these treatments are not mutually exclusive. Hence leading to possible combinations or arms to decide from. We parameterize a total of 28 features for the linear bandits, corresponding to the main effects of the treatments, and the twoway interactions. Given this feature set, each bandit algorithm is seeded with an initial form of experimentation from a run factorial experiment, the minimum balanced design required to estimate all the features.
The response surfaces for experimentation are simulated based off of the general linear model, with main, effects, twoway interactions, threeway interactions, no higher order effects, and noise from the Laplace distribution, as shown in the equation below.^{1}^{1}1Similar experiments were conducted with Gaussian noise. As mentioned in the introduction, the results of the two bootstrap approaches were comparable to the stateoftheart for different noise levels. Due to lack of space we are not reporting the results here.
(1) 
The selection of response surfaces which could contain threeway interactions was intentionally included despite the bandit parameterizations ending with twoway interactions. This was done to reflect the likely scenarios of missing features in any parameterization exercise. One could also think of these as unknown features.
In this study we look at a simulation of response surfaces characterized by different versions of equation 1 that assume different values for the and noise. In order to make these reward generation functions typical of the realworld systems that the bandits are likely to encounter, this study uses a hierarchical probability metamodel (HPM). The HPM, originally proposed in Chipman et al., 18, provides the mathematical structure to determine . This structure seeks to capture certain mathematical regularities seen in realworld response surfaces. Specifically, three properties are sought after. They include, sparsity, hierarchy, and heredity. Sparsity indicates that only a subset of the potential set of features or inputs will be found to have a statistically significant effect on the response. Hierarchy indicates that main effects tend to have larger coefficients than twoway interactions, and twoway interactions tend to be larger than threeway, and so on. Heredity refers to the fact that the likelihood of an interaction having a significant effect on the response is affected by whether it’s parents were significant, or not.
Capturing these properties in the HPM requires modelling various parameters as random variables with conditional dependencies. This can represented as a Bayesian Network (BN)
[Nielsen and Jensen2009]. The structure of using the BN to implement the HPM is illustrated through figure 1. This representation shows the probabilistic relationships in a system with maineffects and their resulting interactions. However, in our simulation we use a system with main effects. The BN helps us determine the statistical significance of the various s, which in turn help us implement the principle of sparsity and heredity. The BN also helps us infer the distributions of the s which helps us implement hierarchy (the coefficient is modelled as a random variable).A meta study by Li et al., 14 adapts this HPM structure in conjunction with realworld data to estimate the parameters shown in Figure 1 . The parameters for the distribution and significance of the coefficients is based off of data sets gathered from published academic work of engineering systems, spanning different disciplines. The adapted mathematical structure of HPMs can be seen in Frey and Li 15. The use of such a test bed to evaluate experimental algorithms can also been seen in other studies [Frey and Li2008, Sudarsanam and Frey2011]. The use of realworld meta data inspires our study to an environment of only arms for the linear bandits problem, as opposed to larger set of arms seen in other studies which use purely synthetic data.
4 Results and Discussion
The pseudoperformance of the selected arm, in line with pseudoregret discussed in [Audibert et al.2009]), is described by:
(2) 
where is the arm selected at trial which is represented as a vector of its feature values, is the optimal arm. The true parameters or feature weights are represented by (which are actually unknown to the agents) and described by their extended feature space of the intercept, the main effects, twoway interactions and threeway interactions ( in total), since they represent the true underlying system. The algorithms (or agent), however, only tries to estimate the partial set of parameters discussed in section 3.
We report our findings across three levels of noise, low (), medium(), and high (), where the noise, is distributed as a Laplace random variable as shown in equation 1.^{2}^{2}2As a reference for why we consider these levels as low, medium and high, the distribution of the significant main effects after scaling of the meta data is and the probability of the main effect being significant is . Hence at the high levels, the strength of the noise is as high as the the significant maineffects, and similar ratios can be derived for medium and low noise levels.
The implementation of XRandom, XFixed, OFUL and LinUCB require making some guesses on tuneable parameters which will influence the performance of the algorithm. The insight on what values are ideal need not be available to the agent, apriori. However, it might be possible to make educated guesses based on the time horizon (number of trials), and degree of noise. To represent these algorithms and their baselines in the best case scenario, their parameters were fine tuned in the simulations. This was achieved with the objective of minimizing cumulative pseudoregret over the 300 trial runs shown in the Figures 2, 3, and 4. In this study, we also wanted to look at the effect of different number of trials, its direct impact on the algorithms as well as its interactive effects with the various noise levels. Table 1 shows this. The cumulative regret shown in the table is an extension of equation 2 defined by
The results show that across all tested degrees of noise and all time horizons, one of the proposed methods, the Random bootstrap, shows lower cumulative regret than all the baseline methods. In certain environments, such as medium to high levels of noise, and lower number of trials, the difference between Random and the comparisons is marginal. However, in other environments, such as low to medium levels of noise, and a larger number of trials, the difference is substantial. The cumulative regret in such cases for the Random is less than of the best baseline approach. In general, the Fixed bootstrap also outperforms the baselines but this is not consistent across all the tested dimensions. Similar to the competitive advantage trend between Random and the baselines, the Fixed also performs relatively better and in lower noise settings and larger number of trials.
In terms of the relative performances between the bootstrap approaches, we see that in general the Random outperforms fixed across the settings that were tested. However, it is noteworthy that Fixed comes closest to the Random for lower trial runs and even outperforms the Random at the lowest noise level for the lower trial setting. As discussed in section 2.2, this is to be expected and is due to the fact that the initial distribution of X is determined from a designed experiment (which is in line with the assumptions of Fixed) and not a random sample. However, subsequent arm selections are the result of inferences from the reward, which is a random draw from a distribution. This favours the XRandom in subsequent trials.
In addition to the comparison of the proposed algorithms with the baselines, the unique testbed proposed in this study could also throw some light on the relative performances of the baseline algorithms when the distributional assumptions used by them are violated. Thompson sampling, when compared to OFUL and LinUCB, consistently performs poorly when the number of trials is small and performs best for longer timehorizons. The poor performance for the shorter trial limits is due to the fixed level of exploration in this algorithm which cannot be finetuned. However, for longer time horizons Thompson sampling outperforms OFUL and LinUCB even after the latter two have had their parameters finetuned. For the greater part, both LinUCB and OFUL perform similarly. This should also be expected owing to the similarities in the algorithms, the distributional assumptions, and the nature of the tuning parameters. However, small differences were seen between the two algorithms, with OFUL working better with fewer trials and lower noise, where as LinUCB did better with higher noise and a longer trial horizon.
The study concludes that when the distribution of noise is not known, and there are no computational constraints, then the Random and fixed bootstraps can be a preferred alternative to currently popular linear bandit algorithms.
4.1 Future Work
Broadly, our future efforts seek to study the role of using distributionfree algorithms to increase the performance and applicability of bandit formulations. To that end, we would like to look at the use of nonparametric statistics techniques to solve various bandit extensions. We would also like to understand the specific situations or environments where the popular methods which make distributional assumptions tend to fail or underperform.
Specifically, our next steps would be to create analytical regret bounds for the linear bootstrap approaches and look at efficient implementations of making the bootstrapped algorithms faster.
References
 [AbbasiYadkori et al.2011] Yassin AbbasiYadkori, David Pal, and Csaba Szepesvari. Improved algorithms for linear stochastic bandits. In Advances in Neural Information Processing Systems, pages 2312–2320, 2011.

[Abe and Long1999]
Naoki Abe and Philip M. Long.
Associative reinforcement learning using linear probabilistic concepts.
In ICML, pages 3–11, June 1999.  [Agrawal and Goyal2012] Shipra Agrawal and Navin Goyal. Thompson sampling for contextual bandits with linear payoffs. arXiv preprint, (1209.3352), 2012.

[Audibert et al.2009]
JeanYves Audibert, Remi Munos, and Csaba Szepesvari.
Explorationexploitation tradeoff using variance estimates in multiarmed bandits.
Theoretical Computer Science, 410(19):1876–1902, Apr 2009. 
[Auer2002]
Peter Auer.
Using confidence bounds for exploitationexploration tradeoffs.
Journal of Machine Learning Research
, 3(1):397–425, 2002. 
[Awerbuch and Kleinberg2004]
Baruch Awerbuch and Robert D. Kleinberg.
Adaptive routing with endtoend feedback: Distributed learning and
geometric approaches.
In
Proceedings of the thirtysixth annual ACM symposium on Theory of computing
, pages 45–53. ACM, June 2004.  [Baransi et al.2014] Akram Baransi, OdalricAmbrym Maillard, and Shie Mannor. Subsampling for multiarmed bandits. In Machine Learning and Knowledge Discovery in Databases, pages 115–131. Springer Berlin Heidelberg, 2014.
 [Breiman and Spector1992] Leo Breiman and Philip Spector. Submodel selection and evaluation in regression. the xrandom case. International Statistical Review/Revue Internationale de Statistique, pages 291–319, 1992.
 [Breiman1992] Leo Breiman. The little bootstrap and other methods for dimensionality selection in regression: Xfixed prediction error. Journal of the American Statistical Association, 87(419):738–754, 1992.
 [Chipman et al.1997] Hugh Chipman, Michael Hamada, and C.F.J. Wu. A bayesian variable selection approach for analyzing designed experiments with complex aliasing. Technometrics, 39(4):372–381, Nov 1997.
 [Dani et al.2008] Varsha Dani, Thomas P. Hayes, and Sham M. Kakade. Visual information extraction with Lixto. In COLT, pages 119–128, July 2008.
 [Eckles and Kaptein2014] Dean Eckles and Maurits Kaptein. Thompson sampling with online bootstrao. arXiv preprint arXiv, 1410(4009), 2014.
 [Efron and Tibshirani1994] Bradley Efron and Robert J. Tibshirani. An Introduction to the Bootstrap. CRC Press, 1994.
 [Frey and Li2008] Daniel D. Frey and Xiang Li. Using hierarchical probability models to evaluate robust parameter design methods. Journal of Quality Technology, 40(1):59, 2008.
 [Kaelbling1994] Leslie Pack Kaelbling. Associative reinforcement learning: Functions inkdnf. Machine Learning, 15(3):279–298, 1994.
 [Li et al.2006] Xiang Li, Nandan Sudarsanam, and Daniel D. Frey. Regularities in data from factorial experiments. Complexity, 11(5):32–45, 2006.
 [Li et al.2010] Lihong Li, Wei Chu, John Langford, and Robert E. Schapire. A contextualbandit approach to personalized news article recommendation. In 19th international conference on World wide web, pages 661–670. ACM, April 2010.
 [Montogomery1984] Douglas C. Montogomery. Design and Analysis of Experiments. Wiley, New York, New York, 1984.
 [Nielsen and Jensen2009] Thomas Dyhre Nielsen and Finn Verner Jensen. Bayesian Networks and Decision Graphs. Springer Science and Business Media, 2009.
 [Osband and Roy2015] Ian Osband and Benjamin Van Roy. Bootstrapped thompson sampling and deep exploration. arXiv preprint arXiv, 1507(00300), July 2015.
 [Rusmevichientong and Tsitsiklis2010] Paat Rusmevichientong and John N. Tsitsiklis. Linearly parameterized bandits. Matehmatics of Operations Research, 35(2):395–411, 2010.
 [Sudarsanam and Frey2011] Nandan Sudarsanam and Daniel D. Frey. Using ensemble techniques to advance adaptive one‐factor‐at‐a‐time experimentation. Quality and Reliability Engineering International, 27(7):947–957, 2011.
 [Sutton and Barto1998] Richard Sutton and Andrew Barto. Reinforcement Learning: An Introduction. MIT Press, Cambridge, Massachusetts, 1998.
 [Thompson1978] Mary L. Thompson. Selection of variables in multiple regression: Part i. a review and evaluation. International Statistical Review/Revue Internationale de Statistique, pages 1–19, 1978.
 [Woodroofe1979] Michael Woodroofe. A onearmed bandit problem with concomitant variable. Journal of American Statistical Association, 74(368):799–806, Dec 1979.