A biobjective method for sample allocation in stratified sampling

A biobjective method for sample allocation in stratified sampling Emilio Carrizosa and Dolores Romero Morales Facultad de Matem´aticas, Universidad de Sevilla, Tarfia s/n, 41012 Sevilla, Spain. email: [email protected]. Sa¨ıd Business School, University of Oxford, Park End Street, Oxford OX1 1HP, United Kingdom; e-mail: [email protected]k. Abstract The two main and contradicting criteria guiding sampling design are accuracy of estimators and sampling costs. In stratified random sampling, the sample size must be allocated to strata in order to optimize both objectives. In this note we address, following a biobjective methodology, this allocation problem. A two-phase method is proposed to describe the set of Pareto-optimal solutions of this nonlinear integer biobjective problem. In the first phase, all supported Pareto-optimal solutions are described via a closed formula, which enables quick computation. Moreover, for the common case in which sampling costs are independent of the strata, all Pareto-optimal solutions are shown to be supported. For more general cost structures, the non-supported Pareto-optimal solutions are found by solving a parametric knapsack problem. Bounds on the criteria can also be imposed, directing the search towards implementable sampling plans. Our method provides a deeper insight into the problem than simply solving a scalarized version, whereas the computational burden is reasonable. Keywords: Integer programming; stratified random sampling; sample allocation; biobjective integer program; parametric knapsack problem. 1 Introduction The sample allocation problem for stratified simple random sampling is the following: we are given a population of size Ndivided into ngroups (strata), with population sizes N1, . . . , Nn. Simple random samples without replacement of sizes x1, . . . , xn,are to be drawn independently from the different strata. The sampling cost within each stratum is assumed to be linear in its sample size xi, with unit sampling cost within stratum iequal to a positive integer ci.The total sampling cost is the sum of the sampling costs within the strata. The drawn sample is used to estimate some parameter of the variable under study Y. Throughout this paper, we assume that the parameter to be estimated is Y , the average of the variable Yin the population. Then, the parameter Ywill be estimated via its Horvitz-Thompson estimator ˆ Y , ˆ Y= n X i=1 Ni Nyi,(1) where yidenotes the sample average within stratum i, see e.g. [5] for further statistical details on the problem considered. Estimator ˆ Yis unbiased, and its variance var(ˆ Y) is given by var(ˆ Y) = Pn i=1 Ni N2var(yi) =Pn i=1 Ni N21 xi−1 Niσ2 c,i,(2) where σ2 c,i is the quasivariance of Ywithin stratum i. We assume, as customary in the literature, that the quasivariances σ2 c,i are either known from previous similar experiments, or replaced by known upper bounds. For instance, if Yi, the values of variable Ywithin stratum i, is a Boolean variable, we can use the upper bound Ni Ni−1 1 4, [5]. The goal is to determine sample sizes x1, . . . , xnminimizing simultaneously •The total sampling cost •The variance of the Horvitz-Thompson estimator ˆ Y . Two types of constraints are imposed. On the one hand, box constraints are considered on the sample sizes xi, li≤xi≤ui,(3) for positive integers li≤ui,for all i= 1,2, . . . , n. Constraints (3) are motivated as follows. First, at least one element must be sampled from each stratum, since, otherwise, the expression (2) is meaningless; moreover, since sampling is without replacement, no more than Niindividuals can be sampled from stratum i. 2 These trivial bounds 1 ≤xi≤Nimay not be sharp enough for practical purposes. Indeed, if we are not only concerned with the variance of the estimator ˆ Y , but also with the variance of the estimators yiwithin the strata, constraints of the form var(yi)≤µi,(4) for µi>0 given, may be imposed. Constraint (4) can also be written as xi≥σ2 c,iNi Niµi+σ2 c,i , which, as asserted, yields a constraint of type (3). On the other hand, the aim of simultaneous minimization of cost and variance may lead to sampling plans in which one of the two objectives attains a low value at the expense of a very high value on the other. To avoid this, we include also in the model target constraints in the form Pn i=1 cixi≤K∗ var(b Y)≤B(5) for positive K∗and B, allowed also to take the value +∞. The problem under consideration is the following biobjective nonlinear integer program min Pn i=1 cixi, var(b Y) s.t. li≤xi≤ui, i = 1,2, . . . , n Pn i=1 cixi≤K∗ var(b Y)≤B xi∈Z, i = 1,2, . . . , n. (6) Define the constants Ai=σ2 c,i Ni N2, i = 1,2, . . . , n B∗=B+Pn i=1 Ai Ni. Then, by (2), (6) yields after erasing additive constant terms min Pn i=1 cixi,Pn i=1 Ai xi s.t. li≤xi≤ui, i = 1,2, . . . , n Pn i=1 cixi≤K∗ Pn i=1 Ai xi≤B∗ xi∈Z, i = 1,2, . . . , n. (PK∗,B∗) In particular, the monotonicity of the criteria implies that, for K∗≥Pn i=1 ciuiand B∗≥Pn i=1 Ai li,constraints (5) are redundant, and (PK∗,B∗) reduces to min Pn i=1 cixi,Pn i=1 Ai xi s.t. li≤xi≤ui, i = 1,2, . . . , n xi∈Z, i = 1,2, . . . , n. (P∞,∞) 3 The set PK∗,B∗of Pareto-optimal solutions of (6), or, equivalently, of (PK∗,B∗), is sought. We recall that a feasible solution x= (x1, . . . , xn) will be Pareto-optimal for (PK∗,B∗) iff no feasible x∗for this problem exists satisfying n X i=1 cix∗ i≤ n X i=1 cixi n X i=1 Ai x∗ i≤ n X i=1 Ai xi , with at least one of the two inequalities above strict. Alternatively, we could be interested in the set of Pareto outcomes, i.e., ( n X i=1 cixi, n X i=1 Ai xi!:x∈ PK∗,B∗). See e.g. [14] for further details on Pareto-optimality in general settings and [7, 16, 18] and the references therein for results and applications to other combinatorial problems. It immediately follows that PK∗,B∗can be obtained from the set P∞,∞of Paretooptimal solutions to (P∞,∞), PK∗,B∗=(x∈ P∞,∞: n X i=1 cixi≤K∗, n X i=1 Ai xi≤B∗).(7) Hence, we can restrict ourselves to the study of (P∞,∞). Although a full description of the Pareto-optimal set for multiobjective integer problems is usually extremely hard, even in the linear case, [7, 16], it turns out that the structure of (P∞,∞) is simple enough to allow us to obtain an easy characterization of P∞,∞,and, by (7), of PK∗,B∗,under certain conditions usually held in practice. Moreover, when such conditions are not fulfilled, standard Branch-and-Bound techniques can be customized to construct PK∗,B∗.As far as the authors are aware, this is the first time this sample allocation problem is directly addressed as a biobjective problem. See Section 3 and [2, 12] for references on related single-objective models. We will illustrate with a real-world database that finding PK∗,B∗provides a deeper insight into the problem than simply solving a scalarized version, whereas the computational burden is reasonable. In what follows we assume that the threshold values K∗and B∗are such that (PK∗,B∗) is feasible. This can be tested by solving, e.g. with the technique described in Section 3.1, the problem min Pn i=1 Ai xi s.t. Pn i=1 cixi≤K∗ li≤xi≤ui, i = 1,2, . . . , n xi∈Z, i = 1,2, . . . , n, and checking whether its optimal value does not exceed B∗. 4 The remainder of the paper is structured as follows. In Section 2 we consider one of the most usual procedures for generating elements of P∞,∞,namely the weighting approach. It turns out that the set of optimal solutions of such problems, the set of supported solutions, can be easily characterized. Sections 3 and 4 address the problem of describing the non-supported Pareto-optimal solutions of (PK∗,B∗).First arbitrary cost structures are considered, and a branch-and-bound algorithm is designed. Finally, the particular case in which the costs are independent of the strata is studied, showing that the supported solutions are the only Pareto-optimal solutions. Numerical experiments with a realworld database are presented in Section 5. The paper ends with a discussion on extensions and lines of further research. 2 Supported solutions for (P∞,∞) A very popular scalarization strategy in Multiple-Objective optimization is the so-called weighting method, in which the objectives are linearly aggregated: ν, 0< ν < 1 is given, and therefore (P∞,∞) is replaced by the scalar problem min (1 −ν)(Pn i=1 cixi) + νPn i=1 Ai xi s.t. li≤xi≤ui, i = 1,2, . . . , n xi∈Z, i = 1,2, . . . , n, or, setting λ:= ν 1−ν∈(0,+∞),by min Pn i=1 cixi+λPn i=1 Ai xi s.t. li≤xi≤ui, i = 1,2, . . . , n xi∈Z, i = 1,2, . . . , n. (8) By varying λin the interval (0,+∞),the set of optimal solutions obtained this way would yield S∞,∞,the so-called set of supported solutions of (P∞,∞). Obtaining the full set (P∞,∞) of supported solutions has important practical consequences. Indeed, if, as frequently done in multiple-objective problems, the final sampling allocation plan is chosen by minimizing a weighted average of the estimator variance and the cost, or, in other words, by solving a problem of type (8), we know that such a plan is not only Pareto-optimal, but also supported. Hence, if, as a preprocessing step, (P∞,∞) is obtained, then an optimal solution to (8) can be obtained from (P∞,∞) by complete enumeration. On the other hand, when the full set of Pareto-optimal solutions is sought, obtaining first the supported solutions and later the remaining Pareto-optimal ones may lead to important savings in computing times, as described for other multiple-objective combinatorial problems e.g. in [17, 18], and shown for this problem in Section 5. Moreover, as described in Section 4, in the important case in which all costs are equal, all Pareto-optimal solutions are supported, thus a description of the supported solutions yields a description of the full set of Pareto-optimal solutions. 5 Two well-known properties of the supported solutions, stated in the following theorem, will be used in the sequel. Theorem 2.1 Let x∗∈ S∞,∞.One has: 1. x∗∈ PK,B for any K≥Pn i=1 cix∗ iand B≥Pn i=1 Ai x∗ i 2. x∗solves the scalar problem min Pn i=1 Ai xi s.t. Pn i=1 cixi≤Pn i=1 cix∗ i li≤xi≤ui, i = 1,2, . . . , n xi∈Z, i = 1,2, . . . , n. Now we address the problem of describing S∞,∞,which, by Theorem 2.1, consists of Pareto-optimal solutions of (P∞,∞).To obtain such description, we first explore the structure of (8) for a choice of λ > 0. As already discussed some years back by Aggarwal [1], (8) is a nonlinear convex separable integer problem with just box constraints, which can be solved analytically. Indeed, consider for each index ithe convex univariate problem min cixi+λAi xi s.t. xi∈R+,(9) where R+denotes the set of non-negative real numbers. Problem (9) has as unique optimal solution xi(λ), xi(λ) = rλAi ci , i = 1,2, . . . , n. Hence, the optimal solutions of min cixi+λAi xi s.t. li≤xi≤ui xi∈Z (10) are given either by the closest feasible point to xi(λ), in case xi(λ) is outside the range [li, ui],or, else, the point(s) in the set {bxi(λ)c,dxi(λ)e},yielding the lowest objective value. For any positive integer k, the objective value in (10) at xi=kis not greater than at xi=k+ 1 iff k(k+ 1) ≥λAi ci.In other words, the set Si(λ) of optimal solutions for (10) will be of the form Si(λ) =            {li},if xi(λ)≤li {ui},if xi(λ)≥ui {bxi(λ)c},if li< xi(λ)< uiand (bxi(λ)c+ 1)bxi(λ)c>λAi ci {dxi(λ)e},if li< xi(λ)< uiand (bxi(λ)c+ 1)bxi(λ)c<λAi ci {bxi(λ)c,dxi(λ)e},if li< xi(λ)< uiand (bxi(λ)c+ 1)bxi(λ)c=λAi ci. (11) 6 Finally, the set S(λ) of optimal solutions to (8) is the Cartesian product of the sets Si(λ) above. Hence, as soon as the parameter λin (8) is provided, a full description of the whole set S(λ) of optimal solutions is at hand. However, in practice, it is not easy to provide a precise value for λ(or ν), [3]. Nevertheless, it is straightforward from the discussion above to obtain a characterization of the set of supported solutions S∞,∞=∪λ>0S(λ).Indeed, one has Theorem 2.2 Given x= (x1, . . . , xn)∈Zn,with li≤xi≤ui∀i, define x= max i:li<xi≤ui xi(xi−1)ci Ai x= min i:li≤xi<ui xi(xi+ 1)ci Ai . Then, x∈ S∞,∞iff x≥x. Proof: By definition, xis supported iff there exists λ > 0 such that xsolves the corresponding problem (8). By (11), for each given i, one has: {λ:xi∈ Si(λ)}=       (0,li(li+1)ci Ai],if xi=li hui(ui−1)ci Ai,+∞,if xi=ui hxi(xi−1)ci Ai,xi(xi+1)ci Aii,if li< xi< ui. (12) Hence, Tn i=1{λ:xi∈ Si(λ)} 6=∅iff x, the highest lower bound in (12), does not exceed x, the smallest upper bound in (12), and the result follows. 2 Expressions (11) and (12) enable us to give an algorithm to describe the set S∞,∞: Theorem 2.3 The set S∞,∞of supported solutions of (P∞,∞)can be written as the union of at most Nsets of the form S1×S2×. . . ×Sn,where each Sjis either a singleton or consists of two consecutive integers, Sj={kj, kj+ 1}. In particular, S∞,∞can be described in O(Nn)time. Proof: By definition, for each λ, the set S(λ) can be written, as the Cartesian product S(λ) = S1(λ)×. . . ×Sn(λ), where each Si(λ) is, by (11), either a singleton or consists of two consecutive integers, and can thus be described in O(n) time. By (12), each set Si(λ) is constant when λvaries in an interval whose endpoints are consecutive elements of the set Ciof critical values of λ, Ci=li(li+ 1)ci Ai ,(li+ 1)(li+ 2)ci Ai ,...,(ui−1)uici Ai.(13) 7 The total number of such critical values for λhas as upper bound Pn i=1(ui−1−li+1) ≤ Pn i=1(Ni−1) = N−n. Hence, O(N) sets of the form S(λ) need to be constructed, yielding an overall time complexity of O(Nn). 2 Remark 2.4 The proof of Theorem 2.3 gives a procedure to describe S∞,∞: Construct, for each i= 1,2, . . . , n, the set Cias defined in (13); for each λ∈Sn i=1 Ci,obtain, following (11), S(λ).Then, Sn i=1 Sλ∈CiS(λ) = S∞,∞.Moreover, without increasing complexity, the output can be obtained sorted by the first or second criterion. Indeed, given 0 < λ1< λ2, and, for j= 1,2,an optimal solution xjfor (8), with λ=λj,it then follows that Pn i=1 Ai x2 i≤Pn i=1 Ai x1 iand Pn i=1 cix1 i≤Pn i=1 cix2 i.Hence, if the list of breakpoints in Sn i=1 Ci is scanned in increasing order of λ, the corresponding set of optimal solutions obtained will appear sorted in nondecreasing and nonincreasing order for the first and second criterion respectively. Since, by construction, each of the nlists Ciin (13) is already sorted, a data structure such as a heap will allow to sort Sn i=1 Ciin O(Nlog n) time. Since the description of S∞,∞given in Theorem 2.3 requires O(Nn),it turns out that a description of S∞,∞with the images sorted by one of the two criteria is obtained in O(Nn) + O(Nlog n) = O(Nn) time. 2 Example 2.5 As a simple illustration, we have considered data on number of employees by area of industrial activity in the region of Andaluc´ıa, Spain, for the year 2000 as available in the database Tempus of the Spanish National Statistics Bureau, INE, [15]. The total number of employees is N= 231,334,grouped into n= 14 industrial activities (strata). The number Niof individuals per stratum is given in the first column of Table 1. The allocation problem is considered under the commonly used assumption that population quasivariances are independent of the strata, and then chosen to be fixed at 1. The lower and upper bounds imposed are the trivial ones: li= 1 and ui=Ni,for all i= 1, . . . , n. We have considered two scenarios for the costs: in the first scenario, we have assumed all the costs cito be equal (and then fixed at 1), and in the latter, different costs ciare randomly associated with the strata, as depicted in the second column of Table 1. By using Theorem 2.3, it turns out that for both cost scenarios S∞,∞is a set with cardinality 231,321,obtained in 0.45 seconds in a AMD Athlon XP 2400+ with 1Gb RAM, running Debian/GNU Linux 3.0 with kernel 2.4.18. The fact that for these two different cost structures the number of supported solutions is the same is due to the arbitrary size of the strata. Figures 1 and 2 depict in the output space the part of S∞,∞in which the first and second objectives do not exceed K∗= 10,000 and B∗= 0.001 for the two above-mentioned cost scenarios. 8 Nici(2nd scenario) 7068 1 53856 1 20450 4 9812 7 10835 9 8331 8 4808 7 20714 8 30142 5 8864 3 9010 4 17519 9 20665 9 9260 9 Table 1: Data of Example 2.5 0.0000 0.0002 0.0004 0.0006 0.0008 0.0010 0.0012 0 2000 4000 6000 8000 10000 12000 costs variance Figure 1: Output of S∞,∞for equal costs 9 Algorithm: Describing PK∗,B∗ Step 0. Set P:= {x∈ S∞,∞:Pn i=1 cixi≤K∗,Pn i=1 Ai xi≤B∗}and set K:= K. Step 1. If K∈nPn i=1 cixi:x∈ S∞,∞,Pn i=1 Ai xi≤B∗o, set K:= K−1 and GoTo Step 4. Step 2. Find one optimal solution xKto (14) with right-hand side coefficient K. Step 3. If xKis the unique optimal solution, then set P=P ∪{xK}.Else, find the set OKof all optimal solutions of (21) with right-hand side coefficient Pn i=1 Ai xK i ,set P:= P ∪OKand K:= Pn i=1 cixK i−1. Step 4. If K≥K, then GoTo Step 2. Else set PK∗,B∗:= Pand STOP. Some technical implementation issues follow. First, by Remark 2.4, one does not need to fully construct in Step 0 the set S∞,∞: if, at some breakpoint λan optimal solution xfor (8) is obtained by (11) with Pn i=1 cixi> K∗(respectively Pn i=1 Ai xi> B∗) then no breakpoint λ0> λ (respectively no λ0< λ) can yield solutions feasible for (PK∗,B∗). For those values of Kfor which the associated problem must be solved in Step 2, the computational burden can be alleviated by doing some simple preprocessing at Step 0 as well as by using some extra information from the cases already studied of this parametric problem. Indeed, (14) is reformulated in (17), as a knapsack problem, to be solved by a branch-and-bound algorithm, in which the LP-relaxation (18) is used as a bounding scheme. The first step performed to solve such LP-relaxation of the root node is to calculate all ratios, ηij ci, and sort them in increasing order. This sorted list of ratios is independent of the parameter K, and can thus be already calculated in Step 0. Moreover, before solving (14), one may already have feasible solutions (obtained from S∞,∞) and upper bounds (the optimal values previously obtained in Step 2 for higher values of K). Step 3 is the hardest part. First, we must check whether the xKobtained in Step 2 is the unique optimal solution to (14). This can be tested in Step 2, by pruning in the branch-and-bound tree only those nodes whose lower bound is strictly worse than the best incumbent. Moreover, we have to keep all optimal solutions of (14) with minimal variance. If instead of finding all optimal solutions to (21) with right-hand side coefficient Pn i=1 Ai xK i ,we just take one of its optimal solutions, the algorithm above will finally stop with a set P ⊆ PK∗,B∗describing the Pareto outcomes. 4 Describing PK∗,B∗.The case of equal costs In this section we address the important particular case of (PK∗,B∗) where the costs ciare independent of the stratum, showing that, contrary to the general case discussed in Section 3, all Pareto optima are supported. By Theorem 2.3, this yields a closed characterization of PK∗,B∗as well as an O(Nn)-time procedure. 16 We have Theorem 4.1 If c1=c2=. . . =cn,then S∞,∞=P∞,∞.Moreover, PK∗,B∗=(x∈ S∞,∞: n X i=1 cixi≤K∗, n X i=1 Ai xi≤B∗).(22) Proof: Any supported solution is Pareto-optimal, so we only need to show the converse. Given x∗∈ P∞,∞, x∗solves (14) for K:= Pn i=1 cix∗ i.The result can then be derived from Theorem 4.1.1. and Section 4.7 of [9]; however, for the sake of self-containedness, a complete proof is derived here. Define, for i= 1,2, . . . , n, and j= 1,2, . . . , ui−li,the Boolean variables y∗ ij as y∗ ij =1,if j≤x∗ i−li 0,else. Then, y∗is optimal for (17), which, since c1=. . . =cn,is equivalent to its continuous relaxation (18). In particular, θ∗≥0 exists such that (θ∗, y∗) is a saddle-point pair for (19) and for its 0 −1 reformulation. Moreover, θ∗>0,since we have strong duality and the optimal value of (14) is strictly positive. By construction of y∗,we have that (θ∗, x∗) is a saddle-point pair for (20), which shows in particular that x∗∈ S(θ∗)⊂ S∞,∞. Then, (22) follows from (7). 2 Under the assumption of equal costs, the set SK∗,B∗of supported solutions of (PK∗,B∗) can also be obtained directly from the set S∞,∞=P∞,∞: Corollary 4.2 If c1=c2=. . . =cn,then SK∗,B∗=(x∈ S∞,∞: n X i=1 cixi≤K∗, n X i=1 Ai xi≤B∗).(23) Proof: We only need to show the inclusion, SK∗,B∗⊆(x∈ S∞,∞: n X i=1 cixi≤K∗, n X i=1 Ai xi≤B∗). Given x∈ SK∗,B∗, x satisfies by construction the constraints, and is Pareto-optimal for (PK∗,B∗). Hence x∈ P∞,∞=S∞,∞,showing the result. 2 17 5 Numerical results In this section we illustrate the two-phase algorithm proposed in Section 3.2 to describe the set of Pareto-optimal solutions for the database presented in Example 2.5. We only consider here the second scenario of costs, since for the first one the set of Pareto-optimal solutions is equal to the set of supported solutions which has been already calculated in Example 2.5. We have imposed budget and accuracy constraints to the sample allocation problem (6). In particular, we have enumerated the set PK∗,B∗, where K∗= 10,000 and B∗= 0.001. In phase 1, 1394 supported solutions were found, which took 0.05 seconds. In phase 2, we found a total number of 3967 non-supported solutions and the computation time was equal to 741.33 seconds. The number of Pareto-optimal solutions is equal to 5361,and the two-phase algorithm took 741.38 seconds. To illustrate the savings reached by applying the knowledge of S∞,∞, we have calculated the same Paretooptimal solutions without performing phase 1 of the algorithm. This took a total of 905.81 seconds, which incurs in a 22.18% increase in the computation time with respect to our two-phase algorithm. We are also interested in extracting some information about the optimization of the linear knapsack problem (17). We may recall that this problem has been solved using a branch-and-bound algorithm. In Figure 4, we have plotted the error bound of the integer solution of (17) available at the root node of the branch-and-bound tree. We observe that this error bound decreases when the right-hand side of the cost constraint, K, increases. Next to the error bound, we have also plotted the computation time per knapsack problem solved, see Figure 5. From this plot we can see that the computation time per knapsack problem solved tends to increase with the parameter K, probably due to the fact that, when the parameter Kincreases, the number of variables with positive value in the LP-relaxation of (17) also increases and the minimum number of nodes which we should inspect to prove optimality in the branch-and-bound scheme of Horowitz and Sanhi increases. Finally, we have observed that the optimal solution of (17) is found in a very early stage of the branch-and-bound tree and most of the nodes are pruned because they are integer or they are not promising. 6 Further research The results obtained in this paper directly apply also to allocation problems for other estimators or other estimation settings, see Carrizosa and Romero Morales [4]. A very challenging extension of the sample allocation model we have dealt with is considering Yas an `-dimensional variable, Y= (Y1, . . . , Y `),thus the `-dimensional parameter Y= (Y1, Y 2, . . . , Y `) must be estimated from a single stratified random sample. The technique developed in this article is still applicable when, for instance, all `variables are 0 −1 or under the very common assumption that the population quasivariance σ2 c,Y j,i of the variables Yjare independent of the strata, see Carrizosa and Romero Morales [4]. The case where none of these reductions can be made will be addressed in a forthcoming 18 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 4000 5000 6000 7000 8000 9000 10000 11000 costs integrality gap (%) Figure 4: Plot of integrality gap for knapsack problem (17) 0.0 0.1 0.2 0.3 0.4 0.5 0.6 4000 5000 6000 7000 8000 9000 10000 11000 costs time to solve knapsack problem (sec.) Figure 5: Plot of time to find non-supported solutions when K∗= 10,000 and B∗= 0.001 19 paper. Acknowledgements This research is supported by grant BFM2002-04525-C02-02 of Ministerio de Ciencia y Tecnolog´ıa, Spain. References [1] Om P. Aggarwal. Bayes and minimax procedures in sampling from finite and infinite populations. The Annals of Mathematical Statistics, 30:206–218, 1959. [2] K.M. Bretthauer, A. Ross, and B. Shetty. Nonlinear integer programming for optimal allocation in stratified sampling. European Journal of Operational Research, 116:667– 680, 1999. [3] E. Carrizosa, E. Conde, F.R. Fern´andez, and J. Puerto. An axiomatic approach to the centdian criterion. Location Science, 3:165–171, 1994. [4] E. Carrizosa and D. Romero Morales. A biobjective method for sample allocation in stratified sampling. METEOR Research Memorandum, RM/03/019, 2003. [5] W.G. Cochran. Sampling techniques. Wiley, New York, 3rd edition, 1977. [6] M. Eben-Chaime. Parametric solution for linear bicriteria knapsack models. Management Science, 42:1565–1575, 1996. [7] M. Ehrgott and X. Gandibleux. An annotated bibliography of multiobjective combinatorial optimization. OR Spektrum, 22:425–460, 2000. [8] D. S. Hochbaum. A nonlinear knapsack problem. Operations Research Letters, 17:103– 110, 1995. [9] T. Ibaraki and N. Katoh. Resource allocation problems. Algorithmic approaches. MIT Press, Cambrigde, Massachusetts, 1988. [10] L. Kish. Optima and proxima in linear sample designs. Journal of the Royal Statistical Society, 139:80–95, 1976. [11] S. Martello and P. Toth. Knapsack problems, algorithms and computer implementations. John Wiley & Sons, New York, 1990. [12] J. Neyman. On the two different aspects of the representative method: the method of stratified sampling and the method of purposive selection. Journal of the Royal Statistical Society, 97:558–606, 1934. 20 [13] C.J. Skinner, D.J. Holmes, and D. Holt. Multiple frame sampling for multivariate stratification. ISI Review, 62:333–347, 1994. [14] R.E. Steuer. Multiple criteria optimization: theory, computation, and application. Wiley, New York, 1986. [15] Tempus database. Instituto Nacional de Estad´ıstica.http://www.ine.es/. [16] E.L. Ulungu and J. Teghem. Multi-objective combinatorial optimization problems: a survey. Journal of Multi-Criteria Decision Analysis, 3:83–104, 1994. [17] E.L. Ulungu and J. Teghem. The two-phases method: An efficient procedure to solve biobjective combinatorial optimization problems. Foundations of Computing and Decision Sciences, 20:149–165, 1995. [18] M. Vis´ee, J. Teghem, M. Pirlot, and E.L. Ulungu. Two-phases method and branch and bound procedures to solve the bi-objective knapsack problem. Journal of Global Optimization, 12:139–155, 1998. 21