title: Discrepancy of stratified samples from partitions of the unit cube
authors: Markus Kiderlen Florian Pausinger
text: 
  
  
  [2010] 11K38, 60C05 (primary), and 05A18, 60D99 (secondary)
  
  Nonparallel Voice Conversion with 
  Augmented Classifier 
  Star Generative Adversarial Networks
  
      Hirokazu Kameoka, 
  
  Takuhiro Kaneko,
  Kou Tanaka,
  and Nobukatsu Hojo
  H. Kameoka, T. Kaneko, K. Tanaka, 
  and N. Hojo are with NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation, Atsugi, Kanagawa, 243-0198 Japan (e-mail: hirokazu.kameoka.uh@hco.ntt.co.jp).
  
  
  This work was supported by 
  JSPS KAKENHI 17H01763 and JST CREST Grant Number JPMJCR19A3, Japan. 
  
      December 30, 2023
  =====================================================================================================================================================================================================================================================================================================================================================================================================
  
  
  
  
  
  
  
  
  
  We extend the notion of jittered sampling to arbitrary partitions and study the discrepancy of the related point sets. Let =(Ω_1,…,Ω_N) be a partition of [0,1]^d and let the ith point in  be chosen uniformly in the ith set of the partition (and stochastically independent of the other points), i=1,…,N. For the study of such sets we introduce the concept of a uniformly distributed triangular array and compare this notion to related notions in the literature.
  We prove that the expected _p-discrepancy, _p(_)^p, of a point set _ generated from any equivolume partition  is always strictly smaller than the expected _p-discrepancy of a set of N uniform random samples for p>1.
  For fixed N we consider classes of stratified samples based on equivolume partitions of the unit cube into convex sets or into sets with a uniform positive lower bound on their reach. It is shown that these classes contain at least one minimizer of the expected    _p-discrepancy. We illustrate our results with explicit constructions for small N.
  In addition, we present a family of partitions that seems to improve the expected discrepancy of Monte Carlo sampling by a factor of 2 for every N.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Setting and main questions
  
  Given a finite set = {_1, …, _N} of N points in [0,1]^d one way to quantify how well-spread these points are, is to calculate the _p-discrepancy 
  
      _p () := ( ∫_[0,1]^d| #(∩[0,  x[)/N - |[0,  x[| |^p  x )^1/p,
  
  of , in which 1≤ p < ∞, 
  #(∩[0,  x[) counts the number of indices 1≤ i ≤ N such that _i ∈ [0,  x[,
  and  |[0,  x[| is the Lebesgue measure of 
  [0, x[ :=∏_k=1^d [0,x_k[
  with x = (x_1, …, x_d); i.e. the _p norm of the so-called discrepancy function.
  For an infinite sequence  the _p-discrepancy _p(_N) is the _p-discrepancy
  of the first N elements, _N, of .
  
  Another important irregularity measure is the star-discrepancy defined as
  
      D^* () :=  sup_ x∈ [0,1]^d | #(∩[0,  x[)/N - |[0,  x[|  |.
  
  The _2-discrepancy is a well studied and understood measure for the irregularities of point sets. We refer to the book <cit.> and the excellent survey <cit.> for further details. In particular, and in contrast to other measures such as the star-discrepancy, it is known how to construct deterministic point sets with the optimal order of magnitude of the _2-discrepancy; see <cit.>. 
  
  For d=2 the optimal order of the _2-discrepancy for finite point sets is known to be 𝒪(√(log N)/N), which already goes back to a result of Davenport <cit.>.
  The optimality of these constructions follows from a seminal result of Roth <cit.> who derived a general lower bound for the _2-discrepancy of arbitrary sets of N points in [0,1]^d; see e.g. <cit.>.
  
  
  
  While deterministic point sets with small discrepancy are widely used in the context of numerical integration, simulations of different real world phenomena may require an element of randomness. The expected discrepancy of a set _N of N i.i.d. uniform random points in [0,1]^d is of order 𝒪(1/√(N)) and as such independent of the dimension.
  For two-dimensional point sets of N i.i.d. uniform random points, we thus also have an expected discrepancy of order 𝒪(1/√(N)) similar to the two-dimensional regular grid (whose discrepancy is known to get worse as the dimension increases).
  
  Randomized quasi-Monte Carlo (RQMC) sampling is a popular method to randomize deterministic point sets; see <cit.> for an excellent introduction. Clever constructions of deterministic point sets, so called quasi-Monte Carlo (QMC) sampling can significantly improve the asymptotic order of integration errors when compared to classical Monte Carlo sampling. RQMC basically takes a deterministic QMC point set as an input and uses a randomisation technique (e.g. a random shift modulo 1 or a so-called digital shift) to generate a new point set, which can be shown to have improved uniform distribution properties compared to Monte Carlo samples, while still enjoying the advantages of being `random' in theoretical analysis; see <cit.>.
  
  
  
  
  
  
  
  The starting point for our work, can also be considered as a basic RQMC technique which was already discussed in <cit.> in a slightly more general form.
  Classical jittered sampling for N=m^d combines the simplicity of grids with uniform random sampling by partitioning [0,1]^d into m^d axis-aligned congruent cubes and placing a random point inside each of them; see Fig. <ref>. Jittered sampling is sometimes referred to as `stratified sampling' in the literature, but we will use the term `stratified sampling' in a more broad sense as outlined below.
  Motivated by recent progress <cit.> the aim of this paper is to take a systematic look at 
  jittered sampling and its extension based on more general partitions =(Ω_1,…,Ω_N) of [0,1]^d. 
  
  We consider stratified sampling, where [0,1]^d is partitioned into N subsets Ω_1,…,Ω_N  and the ith point in  is chosen uniformly in the ith set of the partition (and stochastically independent of the other points), i=1,…,N. 
  If N=m^d and the partition consists of the above mentioned axis-aligned congruent cubes, we obtain jittered sampling as a special case. Besides results for fixed N, we are also interested in the behavior of stratified samples derived from sequences of partitions when N becomes large. 
  
  At this point, we would like to emphasize that sequences of partitions that can be used in stratified sampling are more general than those in Kakutani's splitting procedure and its variants <cit.>. 
  Apart from the obvious difference that these procedures restrict considerations to d=1, the partitions in the present paper need not be nested. This means that the partition in step N+1 is not necessarily obtained as a refinement of the partition in step N; see also the discussion in  Appendix A.
  
  
  It was shown in <cit.> that the asymptotic order of the star-discrepancy of a point set obtained from  jittered sampling is 𝒪(N^-1/2-1/2d). Thus, taking partitions can significantly improve the expected discrepancy of (random) point sets in small dimensions d≥ 2. 
  We are interested in the following main questions:
  
  
    * In which sense are sequences of stratified sample points uniformly distributed as their number N increases? What is the connection to similar notions for partitions in the literature? 	
  
    * Does stratified sampling yield smaller or larger mean discrepancy than Monte Carlo sampling with N 
  i.i.d. uniform random points? Are there  discrepancy notions and assumptions assuring that stratified sampling is strictly better? 	
  
    * Is there a 'best' partition for a given N in terms of a chosen mean discrepancy? 
  
    * Is there a simple family of partitions {^(N)}_N≥ 1 that gives reasonable results for all N and not just for square numbers of points as in the case of classical jittered sampling?
  
    * Can we improve classical jittered sampling with stratified sampling?
  
  Section <ref>  presents our answers to the above together with open questions for future research. In Section <ref> we prove our main theoretical results and illustrate them with examples. Section <ref> introduces and explores an infinite family of partitions and contains more examples as well as numerical results.
  
  
  
  
   §.§ Stratifications and the star discrepancy
   
  By the celebrated result of Heinrich, Novak, Wasilkowski & Wozniakowski <cit.> there exists
  a set of N points in [0,1]^d with
  
      D^*(𝒫) ≤ c √(d/N)     c.
   
  Aistleitner <cit.>, using a result of Gnewuch <cit.>, has shown that one can take c=10. Doerr <cit.> has shown 
  𝔼  D^*(𝒫) ≳√(d/N) for point sets of N i.i.d. uniformly random points indicating that this is the correct order of magnitude. 
  See <cit.> and references therein for the most up-do-date history of improvements on the constant c; the currently smallest value is c=2.4968 derived in <cit.>
  in which it is also shown that (<ref>) with c=3 holds with very high probability when  is a set of N i.i.d. uniformly random points.
  
  Since the best known construction is purely probabilistic, it is natural to ask whether we can improve upon these upper bounds using stratification. Indeed, Aistleitner muses in <cit.> that a thought-out partition may improve the upper bound. 
  Our strong partition principle (Theorem <ref>) shows that the mean _p-discrepancy of stratified sets is strictly smaller than the mean _p-discrepancy of N i.i.d uniform random points and this could lead to a similar result for the star discrepancy using the technique from <cit.>; see also <cit.>. 
  However, the main obstacle in this context is that in order to see the stratification effect in large dimensions, one needs to subdivide the unit cube into (exponentially in d) many subsets and hence one faces a seemingly unavoidable difficulty if one wishes for a result for small N in large dimension. This is also underlined by the discussion on classical jittered sampling in <cit.>, in which it is detailed why jittered sampling gains in effectiveness over purely random points only around N ∼ d^d.
  We believe that stratifications are most useful in small dimensions in which the stratification effect is significant. 
  
  
  
  In which range of d is the stratification effect most significant?
  
  
  We will illustrate the potential of stratifications in the context of star discrepancy with numerical experiments. For fixed (and small) d we expect that it is possible to improve the constant for a uniformly at random scheme with a stratified scheme similar to the case of the _p-discrepancy. For d=2,3,5 we numerically obtain improvements for the families of partitions studied in this paper; see Table <ref>. 
  
  
  
  
  
  
  
  
  
  
  
  § RESULTS
   
  
  
   §.§ Stratified sampling and uniform distribution.
  
  Let d≥ 1 be given. 
  We consider partitions ={Ω_1, …, Ω_N} of the unit cube in ^d into N Lebesgue-measurable sets, i.e.
   
      [0,1]^d=⋃_i=1^N Ω_i,
  
  and the sets do not overlap in the L^1-sense, so Ω_i∩Ω_j is a Lebesgue-null set for all i j in {1,…,N}. It should be emphasized that we prefer this condition to the stronger one requiring pairwise disjoint sets, as we later want to work with closed sets. When the sets Ω_i and Ω_j are convex, then |Ω_i∩Ω_j|=0 is equivalent to saying that Ω_i and Ω_j do not have any interior points in common.
  For the moment, we pose no other geometric conditions on the partition, in particular the sets Ω_i need not be connected.  
   
  
  
  Any partition gives rise to an N-element stratified sample _ of N random points derived from the partition by picking a random uniform (random w.r.t. the normalized Lebesgue measure) point from each Ω_i in a stochastically independent manner. In contrast to stratified sampling in classical sampling theory (see, e.g. <cit.>), we sample only one point in each of the strata.  
   As ground model for comparison we often will consider the set of Monte Carlo  samples _N consisting of N i.i.d. (independent and identically distributed) uniform random points in the unit cube.  
  
  
  For a set ⊂ [0,1]^d of N∈ points in the unit cube let
  
      Z_()=#(𝒫∩ [0,[ )/N
   
  be the proportion of points falling in a test cube [0,[ with ∈ [0,1]^d. 
  For a Monte Carlo sample,  Z_(_N) is a random variable with mean |[0,[|, but Z_(_) need not be unbiased for |[0,[| when  is a partition of the unit cube. 
  We show in Proposition <ref> below  that Z_(_) has mean |[0,[| for all ∈ [0,1]^d
  if and only if the partition is equivolume, that is, if |Ω_1|=⋯=|Ω_N|. 
  This corresponds to the concept of self-weighting stratifications in classical sampling of finite populations: as the samples in all strata are equally large (just one point per stratum), the strata must be equal in size. 
  The assumption of equivolume partitions is often convenient, as it allows us to interpret the mean of _p^p(_) as an integrated centered pth mean; see Equations (<ref>)  and (<ref>), below. Two examples of equivolume partitions for d=2 are illustrated in Fig. <ref>.
  
  
  Now let {^(N)}_N≥ 1 be a sequence of finite partitions of the unit cube and let _^(N)={_1^(N),…,_N^(N)} be the stratified sample associated to the Nth partition. Note that we use capital letters whenever points are random. 
  The fact that partitions for different N need not be related to one another implies that the set of all sampling points forms a triangular array, and we are thus led to define a uniform distribution property for those; see also <cit.>.
  
  
  	A triangular array =(_1^(N),…,_N^(N))_N∈ with points in [0,1]^d is said to be uniformly distributed, if for every cube [ x,  y[ ⊂ [0,1[^d we have
  	
      lim_N →∞#({_1^(N),…,_N^(N)}∩[, [)/N  = |[ x,  y[|.
  
  	
  	
  
  
  A sequence (_i)  in the unit cube is uniformly distributed in the usual sense, if and only if the 
  triangular array (_1,…,_N)_N∈ is uniformly distributed in the sense of Definition <ref>. Hence, this definition generalizes the usual one. As in the classical case, uniform distribution of triangular arrays is equivalent to the weak convergence of the sequence of 'empirical distributions', where the Nth of those distributions sits on the points _1^(N),…,_N^(N) giving equal mass to each of them. In other words,  
  (1/N) ∑_i=1^N f(_i^(N))→∫_[0,1]^df(), as N→∞, for all continuous functions f on the unit cube.  
  
  Proposition <ref> in Appendix A characterizes partitions leading a.s. to uniformly distributed stratified samples using the strong law of large numbers for triangular arrays. 
  	The most important implication of Proposition <ref> is  that stratification sequences of equivolume partitions are uniformly distributed. This is one reason why our theoretical results are based on equivolume partitions. 
  	Appendix A also discusses how Definition <ref> relates to similar concepts in the existing literature.
  
  
  
  
  
  
  
  
  
   §.§ The strong partition principle for stratified sampling.
  
  Discrepancy measures can be used to compare sets of sampling points. In the case of a set  of random sampling points the mean _p-discrepancy _p^p() is often employed, where  denotes the probabilistic expectation. One should correctly call  _p^p() the `mean pth power _p-discrepancy', but we prefer the shorter, slightly misleading form for breviety. 
  
  Certainly, a stratified sample need not be better than a Monte Carlo sample. Consider for instance a partition  with N sets in 
  [0,1]^2 where the N-1 partitioning sets Ω_1,…, Ω_N-1  are all subsets of [δ,1]^2 with some δ∈ ]0,1[. Then the mean _2-discrepancy satisfies
  
      _2(_)^2   ≥∫_[0,δ]^2( 1_[0,[(X^(N)_N)/N-|[0,[|)^2
         =
       δ^4/4N+δ^6/9(1-2/N)≥δ^4/4N>
       _2(_N)^2,
  
  for all δ>(5/9)^1/4≈ 0.86 and N≥2, where the last inequality uses (<ref>). 
  
  In contrast to this, stratified samples from equivolume partitions are never worse than Monte Carlo samples in terms of the mean _2-discrepancy according to the Partition Principle in <cit.>.
  We strengthen this result in two directions showing firstly that stratified samples from equivolume partitions are strictly better, and secondly that _2-discrepancy can be replaced by 
  _p-discrepancy with arbitrary p>1. 
  The main ingredient of our proof is a result by Hoeffding <cit.> stating that 
  among all Poisson-binomial distributions with given mean, the classical binomial distribution is the most spread-out. 
   (Strong Partition Principle) 
  For any equivolume partition  of [0,1]^d with N≥ 2 sets we have 
  
      _p(_)^p< _p(_N)^p
  
  for all p>1. 
  
  The proof of this theorem will be given in  Section <ref>. One can understand (<ref>) as a continuous analog and extension to the statement in finite population sampling theory that self-weighted stratified sampling is always better (in terms of variance) than simple random sampling, both taken with replacement.
  
  
  
  For illustration, we consider the sequence (^(N)_vert) of vertical strip partitions; see  Fig. <ref> (left), but generalized to the d-dimensional case. Direct calculation confirms
  
      _2(_vert)^2 <_2(_N)^2
  
  for all N≥ 2. Both sampling schemes have the same asymptotic order 1/N, but stratified sampling has a better leading constant; see Section <ref> for details. 
  
  
  
  
  
   §.§ Partitions with best average discrepancy.
  
  
  
  For a given N≥ 2, it is an open problem to assure the existence of  a partition 
  whose associated stratified sample has lowest mean discrepancy among all partitions consisting of  N sets. 
  We will show such existence results for certain equivolume partitions. These results are all based on the topological standard argument that a continuous function attains its minimum on a compact set. This requires the choice of a topology on the family  of compact subsets of [0,1]^d. We have chosen the well-established Hausdorff metric, as  is compact in its generated topology. However, both, the extension of this compactness to partitions and the continuity claim of  _p(_)^p as a function of the partitioning sets, require the continuity of the volume functional. We will assure this by assuming certain regularity conditions. More precisely, we assume that there is r>0 such that the sets Ω_1,…,Ω_N⊂ [0,1]^d of the partition have reach at least r, meaning that for any point  with distance less than r from Ω_i there is a unique closest point to  in Ω_i, i=1,…,N; see Fig. <ref>. The class of such sets is very general and contains for instance all compact convex sets in [0,1]^d (as the reach of a closed convex set is infinity). It also contains any given set whose boundary is a piecewise C^2-curve such that its finitely many vertices are 'convex', provided that r>0 is chosen small enough. Let 𝔓_N(r) be the class of all equivolume partitions consisting of sets with reach at least r>0.
  
  Also smaller classes of partitions can be treated. We name here the class 𝔓_N^conv  of equivolume convex partitions, which might be relevant for applications, as all sets constituting a convex partition of [0,1]^d are actually convex polytopes, with the total number of vertices being uniformly bounded when N is given. Hence, convex partitions can be described using finitely many parameters, and thus they are, at least in principle, computationally tractable. 
  The following main result states the existence of equivolume  partitions yielding the best  mean discrepancy Δ from stratification under regularity. Note that the assumptions on the function Δ, which is some given measure of discrepancy, are very weak.
  
  
  Let r>0 and N∈ be given and assume that Δ:([0,1]^d)^N→ is measurable and bounded. Then there exists (at least) one partition ^* ∈𝔓_N(r) such that the corresponding stratified sample _^* minimizes the mean Δ-discrepancy on 𝔓_N(r); i.e.
  		
      min_∈𝔓_N(r)Δ(_) = Δ(_^*).
  
  	A corresponding statement holds true for 𝔓_N^conv.
  
  The standard notions of discrepancy satisfy the assumptions in the above theorem; see the end of Section <ref> for details.
  
  
  	Let r>0, 1≤ p<∞, and N∈ be given. Then there are partitions ^p ∈𝔓_N(r) 
  		and ^* ∈𝔓_N(r) of [0,1]^d such that 
  		
      min_∈𝔓_N(r)_p(_)^p= _p(_^p)^p,
  
  		and 	
  		
      min_∈𝔓_N(r) D^* (_) =  D^*(_^*),
  
  		respectively.
  	
  	Corresponding statements hold true for 𝔓_N^conv.
  
  
  For illustration we will determine the optimal convex partition of [0,1]^2 for N=2 in Subsection <ref>.
  
  Theorem <ref> and its proof  in Section <ref> indicate that the properties of the discrepancy are only of marignal importance while the regularity assumptions on the partitioning sets are crucial in order to show the continuity and compactness. Generally speaking, for such existence statements to hold, we expect that these regularity conditions are not needed, that is, we expect the existence of an equivolume partition of [0,1]^d minimizing a given (rather general) measure of discrepancy. The reasoning for this conjecture is the surmise that minimizers are typically consisting of regular sets. To illustrate this point consider  the simple case d=1, N=2. Clearly, the class of equivolume partitions contains very complicated 
  	pairs of sets, such as fractals, and it is certainly not closed in the Hausdorff-metric nor in the induced _1-metric for indicator functions. But Corollary <ref> in Section <ref>, shows that the unique
  	mean-_2-disrepancy minimizing equivolume partition in this case is simply {[0,1/2],[1/2,1]} 
  	(up to sets of measure zero), which consists of very regular sets.
  
  
  We cannot claim that the equivolume assumption is needed, but it is crucial for our approach, as it avoids that sets in partition sequences shrink to lower-dimensional sets. Existence statements without the equivolume assumption, though very interesting, would require thus substantially different techniques and are beyond the scope of the present paper.
  
  
  
  
  
  
   §.§ Explicit stratification strategies for arbitrary N.
  
  Next, we suggest and motivate a general and versatile construction of partitions for arbitrary N. 
  We define partitions of the unit square generated by parallel lines which are orthogonal to the diagonal of the square. As a special case we consider the partitions _∗^(N) which are equivolume; see Fig.  <ref> (right).
  In Section <ref> we present numerical evidence that stratified samples based on such partitions improve the expected _2-discrepancy of an N-point Monte Carlo sample roughly by a factor of two. 
  As a comparison, we show in Example 1 in Section <ref> that samples based on vertical strip partitions improve an N-point Monte Carlo sample by a factor of 5/3.
  
  Importantly, this construction enables us also to systematically study the role of the equivolume property. 
  In a first step, in Example 3 in Section <ref> we improve the minimal convex equivolume partition for N=2 obtained in Example 2 in Section <ref> by shifting the separating line along the diagonal. In Section <ref> we extend this analysis to the case N=3. We parametrise all such partitions into three sets and determine the minimal partition among them. It turns out that these partitions into three sets have a rich and interesting global structure with respect to their expected discrepancy. 
  
  Finally, we have examples of partitions within this family and for small N that show that it is possible to improve classical jittered sampling by relaxing the equivolume constraint. 
  
  
  
  
   §.§ Conclusions and open questions
  
  In conclusion, our results show that if partitions are needed to generate stratified samples for arbitrary N,
  we suggest to use lines that are orthogonal to the diagonal of the unit square. Within this family it seems that equivolume partitions  _∗^(N) are a reasonably good pick; see Section <ref> for details.
  
  
  Secondly, our examples for N=2,3,4 show that the expected discrepancy can be improved if we drop the equivolume property. This is in line with the results from <cit.> and deserves further attention. It certainly relates to the well-known general observation that the _2-discrepancy exaggerates the importance of points lying close to the origin (see <cit.>).
   
  
  Are there properties of sequences of partitions, other than equivolume, that improve asymptotically the expected discrepancy of Monte Carlo sampling?
  
  
  Our example for N=4 supports the idea brought forward in <cit.> that classical jittered sampling may not give the lowest expected discrepancy for large N.
  
   
  Is there an infinite family of partitions that generates point sets with a smaller expected discrepancy than classical jittered sampling for large N?
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § PROOFS AND EXAMPLES
   
  
  
   §.§ Proofs for Section <ref>
   
  We now give proofs of the results in Subsection <ref> using the notations and notions introduced there. On several occasions we will need the following lemma, which essentially is a reformulation of the fact that a distribution (in the probabilistic sense) is uniquely determined by its cumulative distribution function, also in the multivariate case; see e.g. <cit.>. Indeed, the proof of the following lemma is based on this fact, if the function involved is split into positive and negative part and the total integrals are normalized. 
  
  
  An integrable function  f:[0,1]^d→ is almost everywhere determined if its integrals 
  ∫_[0,]f() are known for almost all ∈ [0,1]^d. 
  
  In other words, ∫_[0,]f()=0 for almost all  ∈ [0,1]^d implies f()=0 for almost all ∈ [0,1]^d. 
  
  We are now in a position to show the announced characterization of equivolume partitions in terms of the unbiasedness of the proportions Z_() in (<ref>). 
   
  
  For  a partition  of [0,1]^d into N Lebesgue sets Ω_1,…, Ω_N of positive volume
  the following three statements are equivalent.
  
  	
    (i)  is equivolume.
  	
    (ii) Z_()=|[0,]| for all ∈ [0,1]^d. 
  	
    (iii) Z_()=|[0,]| for almost all ∈ [0,1]^d.
  
  
  
   
  The bias is 
  
      Z_()-|[0,]|= 1/N∑_i=1^N [|Ω_i∩ [0,]|/|Ω_i|-N|Ω_i∩ [0,]|]= 1/N∑_i=1^N u_i|Ω_i∩ [0,]|,
  
  where 
  
      u_i=1/|Ω_i|-N,     i=1,…,N.
  
  If (i) holds, the vector  =(u_1,…,u_N) is the zero vector and the bias (<ref>) vanishes for all ∈ [0,1]^d. Hence, (i) implies (ii).
  
  Clearly (ii) implies (iii), so it remains to assume (iii) and deduce (i). Assumption (iii) implies
   that (<ref>) vanishes for almost all ∈ [0,1]^d. This implies
  ∫_[0,] f()d=0 for almost all ∈ [0,1]^d, where we have put 
  f=∑_i=1^N u_i 1_Ω_i.
  Lemma <ref> implies f=0 almost everywhere on [0,1]^d, and as  is a partition, =0. Hence the partition is equivolume. 
   
  
  
  
  
  
  
  
   §.§ Proofs for Section <ref>
   
  
  	Let p>1 be given. Using the variable Z_=Z_() from (<ref>) and applying Tonelli's theorem we see that
  	
      _p(_)^p=∫_[0,1]^d(Z_-|[0,]|)^p ,
  
  	so the equivolume assumption and Proposition <ref> yield
  	
      _p(_)^p=∫_[0,1]^d_p(Z_) .
  
  	Here, 
  	
      _p(Y)=|Y- Y|^p
  
  	is the pth centered moment of a random variable Y. The variable  NZ_ is the sum of N independent (but  not identically distributed) Bernoulli variables with success probabilities q_1(),…, q_N(), where 
  	
      q_i()=|Ω_i∩ [0,]|/|Ω_i|=N|Ω_i∩ [0,]|.
  
  	The distribution of NZ_ is usually called Poisson-binomial distribution with N trials and parameter vector ()=(q_1(),…,q_N()). Its mean is ∑_i=1^N q_i()=N|[0,]|. 
  	
  	Setting U_=(#{X_1,…,X_N}∩ [0,])/N with uniform i.i.d. random variables X_1,…,X_N 
  	in [0,1]^d and using similar arguments as above yields correspondingly  
  	
      _p(_N)^p=∫_[0,1]^d_p(U_) d.
  
  	The variable NU_ has a binomial distribution with N trials and success probability |[0,]|. Its mean is therefore coinciding with the mean of NZ_. 
  	
  	We now use the fact that among all Poisson-binomial distributions with given mean, the binomial is the largest one in convex order. This is formalized in <cit.> (see also the paragraph directly after the statement of this theorem) and implies 
  	
      _p(Z_)≤_p(U_)
   with equality if and only if Z_ has a classical binomial distribution, that is, if and only if q_1()=⋯=q_N()=|[0,]|. Integrating (<ref>) with respect to  now yields (<ref>)  if we can exclude the equality case. 
  	
  	
  	
  	Equality in  (<ref>) would imply equality in (<ref>), and thus 
  	N|Ω_i∩ [0,]|=|[0,]|, i∈{1,…,N},
  	for almost all ∈ [0,1]^2. Hence 
  	∫_[0,] 1_Ω_i()d=∫_[0,]1/N d for almost all ∈ [0,1]^d and all i∈{1,…,N}.
  	Lemma <ref> implies 1_Ω_i=1/N,
  	which is not possible as N≥ 2. 
  
  
  It is worth emphasizing the special case p=2 of (<ref>), which has been used more or less explicitly and generally in the existing literature, as it implies that the mean _2-discrepancy can be described as sum of contributions from the individual sample points. Also for p=4 an explicit integral representation can be stated.
  
  	Let  be an equivolume partition and let q_i(), i=1,…,N be defined by (<ref>).  
  Then 
  	
      _2(_)^2=1/N^2∑_i=1^N ∫_[0,1]^d Q_i() ,
  
  with Q_i()=q_i()(1-q_i())^2+q_i()^2(1-q_i())=q_i()(1-q_i()), 
  and 
  
      _4(_)^4=
      1/N^4∑_i=1^N ∫_[0,1]^d R_i()+6/N^4∑_i=1^N∑_j=1  j i^N
      ∫_[0,1]^d Q_i()Q_j(),
  
  where R_i()=q_i()(1-q_i())^4+q_i()^4(1-q_i()). 
  	
  
  
  According to  (<ref>) with p=2, we have 
  	
      _2(_)^2=∫_[0,1]^d(Z_()) ,
  
  	where  Z_() is given in (<ref>). We have already seen that 
  	N Z_()  is the sum of the independent Bernoulli variables Y_i=1_[0,  x[(X_i) with success probabilities q_1(),…, q_N(), respectively, so N^2(Z_()) =∑_i=1^N  q_i()(1-q_i()). This can be inserted into (<ref>) to obtain the first claim. 	
  
  
  	For the second claim, let W_i=Y_i- Y_i, i=1,…,N, and note that 
  	
      N^4_4(Z_())=∑_i_1,…,i_4=0^N   (W_i_1⋯ W_i_4)=∑_i=1^N  W_i^4+4 2∑_i=1^N∑_j=1  j i^N W_i^2 W_j^2.
  
  	As  
  	
      P(W_i=w)={[         q_i(), if  w=1-q_i(),;       1-q_i(),  if  w=-q_i(). ].
  
  	we have W_i^4=R_i() and W_i^2=Q_i(). 
  
  
  
  
  Inserting this into (<ref>) and applying   (<ref>) with p=4 yields the second assertion. 
  
  
  As an application of the previous proposition, we show that an equivolume partition minimizing the  mean _2-discrepancy exists without any further regularity assumptions  when d=1 and N=2. 
  	
  		An equivolume partion =(Ω_1,Ω_2) of the unit interval [0,1] minimizes the  mean _2-discrepancy
  		of its corresponding stratified point set among all equivolume partitions of two sets if and only if Ω_1 coincides up to a set of measure zero with [0,1/2] or 
  		with [1/2,1]. 
  	
  
   
  	Let an equivolume partition =(Ω_1,Ω_2)  of the unit interval be given. It is determined by the measurable set  Ω_1⊂ [0,1] with 1-dimensional Lebesgue measure 1/2, as Ω_2=[0,1]∖Ω_1 (at least up to a set of measure zero). The functions q_i(·) in (<ref>) are thus  q_1(x)=2|Ω_1∩ [0,x]| and q_2=2x-q_1(x) and Proposition <ref> implies 
  		
      _2(_)^2=1/4∫_0^1  2x-4x^2+2q_1(x)(2x-q_1(x)) x=-1/12+
      			1/2∫_0^1  g_x(q_1(x)) x,
  
  		where g_x(q)=q(2x-q). Clearly q↦ g_x(q) is strictly concave for all x∈ [0,1].  
  	
  	
  	
  		It is easy to see that q(x)≤ q_1(x)≤q(x), where 
  		
      q(x)=2|[0,1/2]∩ [0,x]|    and   q(x)=2|[1/2,1]∩ [0,x]|
  
  		for x∈ [0,1]. Hence, there is an α_x∈ [0,1] with 
  		q_1(x)=α_x q(x)+(1-α_x)q(x). Now, (<ref>), the  concavity of g_x and the fact that g_x(q(x))=g_x(q(x))=max{0,2x-1} yield 
  		
      _2(_)^2   ≥ -1/12+
      			1/2∫_0^1  α_x g_x(q(x))+(1-α_x)g_x(q(x)) x
      			
         =-1/12+
      			1/2∫_1/2^1 (2x-1) x=1/24,
  
  		with equality if and only if α_x∈{0,1} holds for almost all x∈ (0,1) due to the strict concavity. As q_1 is continuous, this can only happen when α_x=1 for all x∈ (0,1) or α_x=0 for all x∈ (0,1). These two cases correspond to q_1∈{q, q}, and thus to the two stated choices of Ω_1. 
  	
   
  
  	  
  
  
   §.§ Example 1: Illustration of partition principle
   
  For illustration, we derive the mean _2-discrepancy of an N-point Monte Carlo sample in [0,1]^d. Using 	(<ref>), the i.i.d. property of the sampling points _1,…,_N and 
  
      (1_[0,[(_i))=|[0,[|(1-|[0,[|)
  
  we  obtain 
  
      _2(_N)^2 = ∫_[0,1]^d(1/N∑_i=1^N 1_[0,[(_i) ) 
      =1/N∫_[0,1]^d|[0,[|(1-|[0,[|).
  
  The latter integral equals ∫_[0,1]^d(∏_i=1^d x_i-∏_i=1^d x_i^2) and can be evaluated explicitly. One obtains 
  
      _2(_N)^2 =[1/2^d-1/3^d]1/N.
  
  In particular, for d=2, we get 
  
      _2(_N)^2 =  5/36N.
  
  
  We now compare this with the mean _2-discrepancy of a stratified sample _vert based on the vertical strip partition in Fig. <ref> generalized to arbitrary d≥ 2 by putting
  
      Ω_i={=(x_1,…,x_d)∈ [0,1]^d: i-1/N≤ x_1≤i/N}
  
  for i=1,…,N. The partition =(Ω_1,…,Ω_N) is clearly equivolume. 
    For given =(x_1,…,x_d)∈ ]0,1]^d we let ι̅:=⌊ N x_1 ⌋
  , and obtain for the success probabilities introduced in the proof of Theorem <ref>
  
      q_i=q_i()= N|Ω_i ∩[0,[ | = 
      ∏_j=2^N x_j×
      1     i ≤a̅, 
      
      N x_1 - ι̅     i = ι̅+1,
      
      0     i > ι̅+1.
  
  Due to independence, the relative number of points Z_(_vert) given by (<ref>) has variance
  
      _2( Z_(_vert))    =( Z_(_vert) ) = 1/N^2∑_i=1^N q_i (1-q_i) 
         = 1/N^2( q_ι̅+1(1-q_ι̅+1) + ι̅∏_j=2^N x_j (1-∏_j=2^N x_j) ).
  
  Therefore, (<ref>)  
  yields 
  
      _2(_vert)^2    =1/N^2∫_[0,1]^d[
      N∏_j=1^N x_j-[(Nx_1-⌊ N x_1 ⌋)^2+⌊ N x_1 ⌋]∏_j=2^N x_j^2]  
         = 1/2^dN- 1/3^d-1N^2∫_0^1[(Nx-⌊ N x ⌋)^2+⌊ N x ⌋]  x.
  
  The one-dimensional integral on the right hand side of this chain of equations evaluates to 
  
      ∑_k=0^N-1∫_k/N^k+1/N[(Nx-k)^2+k]  x=1/3+N-1/2=3N-1/6.
  
  Putting things together, we arrive at 
  
      _2(_vert)^2 =[1/2^d-3N-1/2N1/3^d]1/N<
      _2(_N)^2,
  
  where (<ref>) and N≥ 2 was used. This confirms the general result that equivolume stratification is always strictly better than Monte Carlo sampling. It also shows that this stratification scheme has the same asymptotic order (namely 1/N) as Monte Carlo sampling, but a uniformly better leading constant: for instance, when d=2 we get 
   
      _2(_vert)^2=3N+2/36N^2≈3/5_2(_N)^2
  
   for large N. 
  
  
  
  
  
  
  
  
  
  
  
   §.§ Proofs for Section <ref>
   
  
   Fix A⊂^d. We let A and A be the interior and the boundary of A, respectively. For >0 the -parallel set 
  
      A_={x∈^d: inf_y∈ Ax-y≤}
  
  consists of all points with a distance at most  from A.  
  We recall that a set A⊂^d is said to have reach r>0 if for any 0<<r and x∈ A_ there is a point y∈ A such that x-y<x-z for all z∈ A∖{y}. 
  
  The family of non-empty compact sets will be endowed with the Hausdorff metric d_H given by 
  
      d_H(K,K')=inf{>0: K⊂ K_',  K'⊂ K_},
  
  where ∅ K,K'⊂^d are compact. Let  be the family of nonempty compact sets in [0,1]^d, and 
  for r>0 let _r be the subfamily of sets with reach at least r. The latter contains , the family of 
  non-empty compact convex subsets of [0,1]^d. Crucial for our line of arguments is the fact that all three families are compact in the Hausdorff metric. This statement for  is the famous Blaschke selection theorem <cit.>; a proof for  and  _r can be found in <cit.> and 
  <cit.>, respectively. As a reference of convex geometric notions used in this section, we recommend <cit.>.
  
  Importantly, the volume functional is not continuous on  (for instance,  [0,1]^d can be approximated by finite sets in the Hausdorff metric), but it is continuous on both _r and . This can be seen by means of a Steiner-type result stating for K∈_r that
  
      |K_∖ K|=∑_k=0^d-1κ_d-k V_k(K)^d-k
  
  holds for 0≤<r; see <cit.>. Here, κ_d-k is the volume of the Euclidean unit ball in ^d-k, and V_k(K)∈ is the kth total curvature measure (also called intrinsic volume when applied to convex sets). We use repeatedly that K↦ V_k(K) is continuous on _r. 
  These and more results on sets of positive reach can be found in 
  <cit.>; see also the survey <cit.>, where an outline of the history, newer results and 
  additional references on the matter can be found. 
  
   
  	Let N≥1 and r>0 be fixed. 	
  	The family 𝔓_N(r) of all equivolume  partitions of [0,1]^d consisting of N sets in _r is compact. 
  	
  	The same holds true for the family 𝔓_N^conv of all equivolume partitions of [0,1]^d consisting of N convex sets. 
  
  
  	 Clearly, the family of equivolume partitions of sets in _r  is a subset of the Cartesian product _r^N, 
  	
  	
  	
  	
  
  
  
  more precisely, 
  
      𝔓_N(r)={(K_1,…,K_N)∈_r^N: ⋃_i=1^K K_i=[0,1]^d,  and  
      |K_1|=⋯=|K_N|=1/N}.
  
  In fact, assume that (K_1,…,K_N) is an element of the right hand side of (<ref>). If there was a set K_j overlapping with ⋃_i j^N K_i, we would have 
  
      0<|K_j∩⋃_i j^N K_i|= 
         |K_j|+ |⋃_i j^N K_i|-|K_j∪⋃_i j^N K_i|≤1/N+(N-1)1/N-1=0,
  
  a contradiction.
  
   
  
  From the definitions it is clear that (K,M)↦ K∪ M is Lipschitz continuous if the product space is endowed with the maximum norm of the marginal metrics. 
  Hence, the right hand set in (<ref>) is closed in _r^N, as it is defined by means of continuous functions. As _r^N is compact, 𝔓_N(r) is compact too. 
  
  Finally, ^N is compact implying the compactness of 𝔓_N^conv=𝔓_N(r)∩^N. 
  
  
  
  
  
  
  We now show that an average discrepancy of a stratified sample is continuous as a function of the partitioning sets. 
  In the following, Δ stands for a measure of discrepancy, and can be the _p-discrepancy or any other, which satisfies the rather weak assumptions below. 
  
  
  	
  	Fix r>0. If  Δ:([0,1]^d)^N→ is measurable and bounded, 
  	then ϕ_Δ:_r^N→ with 
  	
      ϕ_Δ(K_1,…,K_N)=∫_([0,1]^d)^NΔ(_1,…,_N)
              ∏_i=1^N 
              _K_i(_i)
           (_1…,_N)
  
  	is Lipschitz continuous. 
  
  
  	We let ·_∞ be the L^∞-norm of bounded functions  on ([0,1]^d)^N. We will use the bound 
  	
      |∏_i=1^N t_i-∏_i=1^N t_i'|≤∑_j=1^N |t_j-t_j'|,
  
  	which holds for all t_1,…,t_N,t_1',…,t_N'∈{0,1}.  Let 0<<1. 
  	  If  d_H(K_i,K_i')≤ for all i=1,…,N, the indicators _K_i and _K_i' coincide on 
  	  the complement of 
  	
      M_i=
      	[(K_i)_∖K_i]∪[(K_i')_∖K_i'].
  
  	By Steiner's formula (<ref>), we have 
  	
      |M_i|≤∑_k=0^d-1^d-kκ_d-k [V_k(K_i)+V_k(K_i')].
  
  	As V_k is continuous and _r is compact, V_k(K_i)+V_k(K_i')≤ 2max_M∈_r V_k(M)<∞, so 
  	
      |M_i|≤ c
  
  	for all i=1,…,N, 
  	with a constant c that does not depend on (K_1,…,K_N,K_1',…,K_N').  
  	Thus, by (<ref>) and (<ref>), we have 
  	
      |ϕ_Δ(K_1,…,K_N)-ϕ_Δ(K_1',…,K_N')|
      	
      	≤ cN Δ_∞.
  
  	This implies the claimed continuity.
  
  
  
  
  
  	Let Δ be as in Proposition <ref> and fix r>0. For any finite equivolume partition =(Ω_1,…,Ω_N) of [0,1]^d with sets in _r, 
  	let _={_1,…,_N} be the corresponding stratified sample. 
  	
  	Then Δ(_1,…,_N) is continuous as a function of ∈𝔓_N(r). 
  
  
  	As the partitions are equivolume, we have |Ω_i|=1/N for all i, so 
   
      Δ(_1,…,_N)
       =N^Nϕ_Δ(Ω_1,…,Ω_N),
  	
  and the claim follows from Proposition <ref>. 
  
  
  
  	Assume that Δ:([0,1]^d)^N→ is measurable and bounded. Proposition <ref> thus implies that Δ(_) is a continuous function of ∈𝔓_N(r), where   _ is the corresponding stratified sample.
  	
  	As 𝔓_N(r)  and its subset 𝔓_N^conv  are both compact by Proposition <ref>,  Δ(_) attains minima on either set. This shows the assertion. 
  
  
  Corollary <ref> now follows directly from Theorem <ref>. In fact, for 1≤ p<∞ the function 	
      Δ_p(_1,…,_N)=_p({_1,…,_N})^p
   is bounded.  
  	It is also continuous due to a dominated convergence argument. Even simpler, the measurability and boundedness of Δ(_1,…,_N)=D^*({_1,…,_N}) follows directly from the definition. 
  
  
  
  It should be remarked that the above arguments do not rely on the choice of the unit cube as reference set. 
  The results and proofs continue to hold with minor modifications if the  set [0,1]^d is replaced by any fixed compact convex set K⊂^d.  
  
  
  
  
  
  
  
  
  
  
   §.§ Example 2: Convex equivolume partitions into two sets
  
  Recall that 𝔓_N^conv is the family of all convex equivolume partitions of [0,1]^d with N elements. According to Theorem <ref>, there exists a partition in 𝔓_N^conv  that minimizes the mean _2-discrepancy. The following result determines this partition for N=2 and d=2. In this simple case, 𝔓_N^conv  can be described with the help of a one-parameter model, which is relatively easy to analyze. The optimal partition _*^(2) is obtained by cutting [0,1]^2 into two congruent triangles by the anti-diagonal; see Fig. <ref>, right. 
  
  
  
  
  
  
  
  
   
  For d=2 we have 
  
      min_∈𝔓_2^conv_2^2(_) = 
      	_2^2(__*^(2) ) = 0.05.
  
  
  
  
  Let Ω_1 and Ω_2 be two convex sets  that partition [0,1]^2 and have the same content. 
  By convexity, the intersection of  Ω_1 and Ω_2 is contained in a line ℓ.  
  The midpoint p=(1/2,1/2) of [0,1]^2  must be contained in one of these sets, so let us assume p∈Ω_1. 
  As the reflection ℓ' of ℓ at p is parallel to ℓ, the reflection Ω_1' of Ω_1 at p must contain Ω_2. By the equivolume property we have |Ω_2|=|Ω_1|=|Ω_1'|, so Ω_2=Ω_1' up to a Lebesgue-null set. As Ω_2 and Ω_1' are closed and convex we must even have Ω_2=Ω_1', so 
  Ω_2 is the reflection of Ω_1 at p and we conclude p∈ℓ. 
  
  As  _*^(2) and _2^2(_) are unaltered if all sets in the partition are reflected at the main diagonal  of [0,1]^2, we may assume from now on that ℓ hits the x-axis in a point A∈ [0,1]; see Fig. <ref> (left). Note that _*^(2) corresponds to A=1. For fixed A, we assume from now on that  Ω_1 is the partitioning set that contains the left, vertical edge of the unit square.
  
  To calculate the expected _2-discrepancy, we use the formula 
  
      _2^2 (_ )    = 1/72 + 2 ∫_[0,1]^2f() x_1x_2 - f()^2  
         = 1/72 + 2 ∫_Ω_2f() x_1x_2 - f()^2
  
  from <cit.>, in which f()= | Ω_1 ∩ [0,] |, and where we used for the second equality that the integrand vanishes whenever ∈Ω_1.
  
  We distinguish the three cases A=[0,1/2), A=(1/2,1] and A=1/2.
  The special case A=1/2 gives the vertical strip partition for N=2 with expected discrepancy 1/18=0.055… according to (<ref>).
  
  Now assume that A=(1/2,1]. In this case, the separating line  ℓ has the equation
  
      y = 1/1-2A x - A/1-2A.
  
  We partition Ω_2 into the two sets Ω_2^+={(x,y)∈ [0,1]^2: x≥ A} and Ω_2^-=Ω_2∖Ω_2^+; see Fig. <ref> (middle). For ∈Ω_2^+  we have
  
      f(x,y) 
      =Ay-(2A-1)y^2/2,
  
  so
  
      2 ∫_Ω_2^+f() x_1x_2 - f()^2   
       =1/120 (-2A^3-11A^2+10A+3),
  
  whereas for ∈Ω_2^- we have
  
      f(x,y) 
      
      =(1-2A)y^2/2+Ay+(A-x)^2/2(1-2A),
  
  which results in
  
      2 ∫_Ω_2^-f() x_1x_2 - f()^2   
      	=1/360 (1-2A)^2  (11+2A).
  
  In total, we obtain for A ∈ (1/2,1]
  by inserting the contributions of both cases into (<ref>) that 
  
      _2^2(_ ) = 1/360(2A^3 +3A^2-12 A+25),
  
  which attains  its minimum for A=1 with a value of 1/20=0.05; see Fig. <ref>.
  Note that as A→ 1/2 this function approaches 0.055….
  
  
  We can analyse the last case A=[0,1/2) in a similar fashion and obtain 
  
      _2^2(_ ) = 1/360(-6A^3+3A^2-6A+23),
  
  which approaches its infimum as A→ 1/2 with a value of 0.055…; see Fig. <ref>. This proves the assertion.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § AN INFINITE FAMILY AND NUMERICAL RESULTS
   
  
  
  
  
   §.§ An infinite family and special cases
  
  Motivated by the result of the previous section we define an (N-1)-parameter family of partitions of the unit square generated by parallel lines which are orthogonal to the diagonal of the square. For N∈ and a vector v=(v_1,…,v_N-1)∈ [0,√(2)]^N-1 with 0<v_1<v_2<⋯<v_N-1<√(2) we define a partition _𝐯^(N) as follows. If ℓ_i denotes the line with slope -1 hitting the first closed quadrant and with distance v_i from the origin, then [0,1]^2∖{ℓ_1,…, ℓ_N-1} has N connected components. Its closures are denoted by Ω_1,…Ω_N, where Ω_i is positioned between ℓ_i-1 and ℓ_i if we use the convention v_0=0 and v_N=√(2).  Examples are illustrated in Fig. <ref>.
  We will often use the abbreviation _v_1,…,v_N-1^(N) for the partition _(v_1,…,v_N-1)^(N). 
  
  
  An interesting special case are the equivolume partitions denoted by _∗^(N). They are defined via
  _∗^(N)= _v_1,…,v_N-1^(N) with
  
      v_i =  √(i/N),
  
  for 1 ≤ i ≤⌊ N/2 ⌋ and 
  
      v_i = √(2) - √(N-i/N),
  
  for ⌊ N/2 ⌋ +1 ≤ i ≤ N-1.
  
  As a side remark, we mention also the partition generated by a set of equidistant points; i.e. v_i = √(2) i/N for a given N>1. This is a simple example of a family of partitions that is not equivolume for any N. However, by pairing complementary sets such that the union has volume 2/N, it is possible to obtain equivolume partitions for every even N with N/2 points into non-connected sets.
  
  
  
  
   §.§ Example 3: Relaxing the volume constraint I
   
  We have seen in Example 2 that _∗^(2) gives the lowest expected discrepancy among all convex equivolume partitions into two sets. We will now show that this partition can be improved if the equivolume condition is dropped by shifting the line along the diagonal, that is, by considering partitions in the class _v^(2) with v∈ [0,√(2) ]; see  Fig. <ref>.
  
  
  We have
  
      min_v ∈ [0,√(2)]_2^2(__v^(2) ) = _2^2(__v^∗^(2) ) ≈  0.049,
  
  for v^∗= 0.793398…. 
  
  A comparison with Lemma <ref> shows that _v^∗^(2) has a smaller mean _2-discrepancy than any convex equivolume partion when N=2.
  
  For a given point v∈ [0,√(2)] we denote the corresponding intersection of the line with the boundary of the square with (A,1) if v ∈ [√(2)/2,√(2)] and with (B,0) for v ∈ [0,√(2)/2].
  If v ∈ [√(2)/2,√(2)], then A=√(2)v-1. The two partitioning sets have volumes 1-(1-A)^2/2 and (1-A)^2/2 and points on the separating line satisfy y=-x+1+A. In total, we obtain for A ∈ [0,1] and in a similar fashion as in the proof of Lemma <ref> that
  
      _2^2 (_A) = -18 - 30 A + A^2 - 36 A^3 + 52 A^4 - 12 A^5 - 2 A^6/-360 - 720 A + 360 A^2.
  
  This function attains its minimum for A=0.122034 with value 0.04904…; see Fig. <ref>. 
  
  Next, if v ∈ [0,√(2)/2], then B=√(2)v. The two partitioning sets have volumes B^2/2 and 1-B^2/2 and points on the separating line satisfy y=-x+B. Similar to the previous considerations we split the integral into subintegrals and distinguish the different cases on which the intersections can be described with the same function. We obtain
  
      _2^2 (_B) =
       -135 + 120 B + 175 B^2 - 288 B^3 + 112 B^4 - 2 B^6/360 (B^2-2)
  
  for B ∈ [0,1], which attains its minimum for B=1 with value 0.05; i.e. the minimum is attained for the anti-diagonal; see Fig. <ref>.
  This concludes the proof.
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Example 4: Relaxing the volume constraint II
   
  In this section we extend the results of Examples 2 and 3 to the case N=3. The partitions can still be analysed explicitly and it turns out that there is a unique partition that minimises the expected discrepancy; see Fig. <ref>.
  However, the full analysis consists of a tedious case-by-case study and, hence, we do not fully outline the proof of this assertion in the following, but report the most interesting cases only.  We leave the analysis of the remaining cases – which we carried out and which follows along the same lines as our proof – to the interested reader.
  
  In general, let 0< v_1 < v_2 < √(2) and denote the three sets of the partition with Ω_1, Ω_2 and Ω_3, as described in Subsection <ref>. Associated to the three sets, there are three indicator functions χ_1, χ_2 and χ_3 with
  
      χ_j(x,y) := 
      
      1     if X_j ∈ [0,x]× [0,y]
      , 
      
      0     else,
  
  where X_j is the random point in Ω_j. 
  
  Setting 
      #(x,y) := χ_1(x,y) + χ_2(x,y) + χ_3(x,y),
  
  we get for the expected _2-discrepancy
  
      ( #(x,y)/3 - x y  )^2 = 1/9(#(x,y) )^2 - 2/3 x y (#(x,y)) + x^2 y^2,
  
  and since #(x,y) is a Poisson-binomial distributed random variable we have that
  
      (#(x,y)) = q_1(x,y) + q_2(x,y) + q_3(x,y),
  
  in which, as in (<ref>),
  
      q_j(x,y) = (χ_j=1) =  | [0,x]× [0,y] ∩Ω_j |/|Ω_j|.
  
  The analysis proceeds now via two levels of case distinctions. 
  The first level concerns the actual partitions into three sets that we are considering; see Fig. <ref>.
  Within each of the four cases, we partition the unit square in dependence of v_1 and v_2 into sets in which the probabilities q_j have the same closed form in x,y, thus providing closed formulas for (#(x,y)) and (#(x,y))^2; see Fig. <ref>.
  Next, we compute the expected discrepancy, which is an integral over the unit square, as a sum of integrals over the different closed form expressions according to the subcase-partition.
  
  
  
  
  
  
  
  In the following, we first focus on the third case in Fig. <ref>. In this case we have  0 < v_1 < 1/√(2) < v_2 < √(2) with v_2 ∈ [2v_1, √(2)] and A= v_1√(2), B= v_2√(2)-1 such that for given A we have that 2A-1 ≤ B ≤ 1. 
  It turns out that we can calculate the expected discrepancy as a rational function of A and B, which has a unique minimum.
  For simplicity we restrict considerations to 1/2≤ A ≤ 1 
  
  in Lemma <ref>.
  
   
  Let 1/(2 √(2)) ≤ v_1 < 1/√(2) and 2v_1 ≤ v_2 ≤√(2)-√(2)-2v_1/2.
  
  Let _v_1,v_2^(3) be the corresponding partition of the unit square into three sets. Then
  
  
      min_v_1,v_2_2^2(__v_1,v_2^(3) ) = _2^2(__v_1^∗,v_2^∗^(3) ) =0.0267804…,
  
  for v_1^∗ = 0.5130… and v_2^∗ = 1.1249…
  
  
  
  We set A= v_1√(2), B= v_2√(2)-1 such that for given 1/2 ≤ A ≤ 1 we have 2A-1 ≤ B ≤ A.
  
  We get |Ω_1|= A^2/2, |Ω_3| = (1-B)^2/2 and |Ω_2| = 1-|Ω_1| - |Ω_3|. 
  To calculate the expected discrepancy we subdivide the unit square into 6 sets, S_1, …, S_6, as illustrated in Fig. <ref> and we use the symmetry along the diagonal in cases II-VI.
  For a point (x,y) ∈ S_i we denote the expected discrepancy function by f_i(x,y)   and we calculate the expected discrepancy for the whole partition as
  
      _2^2(__v_1,v_2^(3) ) = ∫_S_1 f_1(x,y)  (x,y) + 2 ∑_i=2^6 ∫_S_i f_i(x,y)  (x,y).
   
  
  For illustration, we give the calculations for the first case and refer to Appendix B for the other (equally elementary, but much more technical) cases.
  
  
    
  Case I Let (x,y) ∈Ω_1=S_1. Then q_2(x,y)=q_3(x,y)=0 and q_1(x,y)=2xy/A^2. Hence, we get
  
      f_1= ( #(x,y)/3 - x y  )^2 = x y(2 + 3 (-4 + 3A^2)x y)/9 A^2,
  
  and
  
      ∫_0^A∫_0^A-xx y(2 + 3 (-4 + 3A^2)x y)/9 A^2 y  x = 1/540 A^2 (5-4A^2 + 3A^4).
  
  
  
    
  Cases II - VI With similar calculations as in Case I we can obtain expressions f_2, …, f_6 for the expected value of the discrepancy function for points (x,y) in each of the sets and in dependence of the parameters A and B; see  Appendix B for details. Integrating these functions over their respective domains and summing the values, 
  
  
  
  
  
  gives us the following rational function in A and B which we can minimize over 1/2≤ A ≤ 1 and 2A-1 ≤ B ≤ A in order to obtain the two parameters A and B that generate the partition with the smallest expected discrepancy in this family; i.e. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
      _2^2(__v_1,v_2^(3) ) =
       1/12960 A^2 (-1 + B)^2 (-1 + A^2 + (-2 + B) B)×
          (128 B^7 - 64 B^8 + A^10 (-1440 + 2880 B - 1440 B^2)  . 
         +A^9 (2592 - 2592 B - 2592 B^2 + 2592 B^3)  
         + A^8 (-3832 + 3936 B - 1248 B^2 + 5760 B^3 - 4680 B^4)  
         + A^7 (10368 - 12288 B - 4032 B^2 + 8640 B^3 - 6336 B^4 + 4032 B^5)  
         + A^6 (-13936 + 19408 B + 4360 B^2 - 12864 B^3 - 1104 B^4 + 6480 B^5 - 
          3240 B^6)  
         +A^5 (5472 - 8544 B - 8736 B^2 + 20384 B^3 - 5280 B^4 - 2400 B^5 - 
          1440 B^6 + 1440 B^7)  
          +A^4 (1908 - 1920 B + 5496 B^2 - 8976 B^3 - 480 B^4 + 5760 B^5 - 
          1608 B^6 - 144 B^7 - 36 B^8) 
            +  A^3 (-480 - 1440 B + 1920 B^2 + 3840 B^3 - 7840 B^4 + 3104 B^5)  
          +A^2 (-787 + 1740 B + 528 B^2 - 6292 B^3 + 9198 B^4 - 3492 B^5 
            -212 B^6 + 84 B^7 + 237 B^8 - 72 B^9 - 36 B^10)
         .  + A (-768 B^6 + 384 B^7) )
  
  This function can be minimised using a standard computer algebra system.
  The minimum of this function is 0.0267804 for the parameter values A= 0.725501 and B= 0.590843. 
  These parameter values correspond to v_1 = 0.5130… and v_2 = 1.1249….
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  In a similar fashion we can analyse the second case in Fig. <ref>.
  
  
   
  Let 1/(2 √(2)) ≤ v_1 < 1/√(2) and 1/√(2)≤ v_2 ≤ 2 v_1.
  
  Let _v_1,v_2^(3) be the corresponding partition of the unit square into three sets. Then
  
  
      min_v_1,v_2_2^2(__v_1,v_2^(3) ) = _2^2(__v_1^∗,v_2^∗^(3) ) =0.0268054…,
  
  for v_1^∗ = 0.5329… and v_2^∗ = 1.06582…
  
  
  
  The proof follows the same lines as in Lemma <ref>; the only difference is the subdivision of the unit square as well as the range of B.
  We set again A=v_1 √(2), B=v_2 √(2)-1 such that for given 1/2 ≤ A ≤ 1 we have 0 ≤ B ≤ 2A-1.
  
  We get |Ω_1|= A^2/2, |Ω_3| = (1-B)^2/2 and |Ω_2| = 1-|Ω_1| - |Ω_3|. 
  To calculate the expected discrepancy we subdivide the unit square into 6 sets, S_I, …, S_VI, as illustrated in Fig. <ref> and we use again the symmetry along the diagonal in cases II-VI.
  For a point (x,y) ∈ S_i we denote the expected discrepancy function  by f_i(x,y) and we calculate the expected discrepancy for the whole partition as
  
      _2^2(__v_1,v_2^(3) ) = ∫_S_1 f_1(x,y)  (x,y)  + 2 ∑_i=2^6 ∫_S_i f_i(x,y) (x,y).
  
  
  
  
  
  
  
  As before, explicit expressions for f_i on S_i, i=1,…,6 can be obtained. They depend on the parameters A and B. Integrating these functions over their respective domains and summing the values, 
  
  
  
  
  
  gives us again a rational function in A and B which we can minimize over 1/2≤ A ≤ 1 and 0 ≤ B ≤ 2A-1 in order to obtain the two parameters A and B that generate the partition with the smallest expected discrepancy in this family. More explicitly, we have
  
      _2^2(__v_1,v_2^(3) ) =
       1/1620 A^2 (B-1)^2 (A^2+(B-2) B-1)×
           ( A^8 (-6 B^2+12 B-14)+48 A^7 B+A^6 (-3 B^4+12 B^3+140 B^2-544 B+283) . 
            +A^5 (208 B^3-672 B^2+1152 B-576) 
            + A^4 (-3 B^6-18 B^5-15 B^4-420 B^3+1065 B^2-942 B+333) 
            +A^3 (208 B^5-440 B^4+480 B^3-480  B^2+120) 
            +A^2 (-6 B^8-24 B^7+134 B^6-444 B^5+870 B^4-704 B^3+264 B^2+168 B-146) 
           . +A (48 B^7-96 B^6)-8 B^8+16 B^7 ).
  
  Using again a computer algebra system we obtain that the minimum of this function is 0.0268054 for the parameter values A= 0.753647 and B= 0.507294. These parameter values correspond to v_1 = 0.5329… and v_2 = 1.06582….
  
  
  
  
  
  
  
  
  
  
  The two lemmas provide two interesting insights. 
  Firstly, we can now easily analyse the equivolume partition within this family.
  
  
  Let v_1 = √(1/3), and v_2 = √(2) - √(1/3), then 
  the mean _2-discrepancy of _∗^(3)= _v_1,v_2^(3) is 
  
      _2^2(_Ω_∗^(3) ) = 0.0290077.
  
  
  
  We have that A= √(2)/√(3) =0.816… and B= 1 - √(2)/√(3) = 0.183….
  Hence, we see that this case satisfies the assumptions of Lemma <ref>. Using (<ref>)  we obtain the value.
  
  
  Secondly, combining the results of the two lemmas, we can now fix a parameter A in [1/2,1] and analyse all partitions for this fixed A and any parameter B in [0,A]. 
  As it turns out, if we fix A and plot the expected discrepancy as a function of the parameter B, then this function is very well behaved and has one unique minimum; see Fig. <ref>. 
  
   
  It is interesting to note that the minimal parameters in Lemma <ref> are not at the boundary of the two cases; i.e. for A= 0.72550 we have that 2A-1 = 0.451003 < 0.590843 = B < 0.72550 = A. The expected discrepancy for the partition generated by (A,2A-1) is 
  
      0.0269763
  
  and is thus only slightly larger.
  Furthermore, it turns out that the minimal parameters in Lemma <ref> are exactly at the boundary; i.e. for A= 0.753647, we have that 2A-1=B= 0.507294, whereas the minimum for this A is obtained for B^*=0.516474 and the expected discrepancy is
  
      0.0268046
  
  and is thus only slightly smaller.
   
  
  Interestingly, the minimum for a given A can be obtained in either of the two cases analysed in Lemma <ref> and <ref> as can be seen from the examples. As a rule of thumb, the global minimum within this family is obtained for parameters A and B in which the minimum for fixed A lies almost at the interval boundary, i.e. for which B_min≈ 2A-1.
  
  This observation relates to Question <ref> of Section <ref>. It illustrates that the equivolume property appears to have no particular significance within this simple family of partitions. It rather seems that other geometric reasons drive the minimisation.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ An algorithmic approach
  
  In order to run systematic experiments within this family of partitions, we implemented an algorithm that takes as input an arbitrary vector 𝐯={v_1, …, v_N-1} with increasing entries in [0,√(2)] as well as a point (x,y) ∈ [0,1]^2 and outputs the expected value of the discrepancy function of the set of N points generated from the partition _𝐯^(N) on the interval [0,x] × [0,y]. This allows for an approximation of the expected value of the _2-discrepancy of __𝐯^(N) using standard results from the theory of quasi-Monte Carlo integration. 
  
  The algorithm is based on a simple geometric consideration. Assume 0 ≤ y≤ x ≤ 1. As we have seen, 
  we need to determine the probability 
  
      q_i=q_i(x,y)= P(_i ∈Ω_i ∩[0,x] × [0,y]) = | Ω_i ∩ [0,x] × [0,y] | /| Ω_i |
  
  with which the point _i ∈Ω_i lies in the box [0,x] × [0,y] for 1≤ i ≤ N. The expectation is then obtained from Proposition <ref>. 
  Hence, we first need to calculate the respective areas of the sets Ω_i. To calculate their intersections with [0,x] × [0,y] we divide the set v_1,…,v_N-1 into four subsets depending on which of the vertices of [0,x] × [0,y] are on the left or on the right of ℓ_1,…,ℓ_N-1, respectively. More precisely:  the four lines with slope -1 through the vertices of [0,x] × [0,y] have distances 0=u_0≤ u_1≤ u_2≤ u_3≤√(2) from the origin. The jth subset consists then of all v_i's between u_j-1 and u_j for j=1,…,4, where  we have put u_4=√(2); see Fig. <ref> (left).  Different formulae are used to compute the intersection in each case.
  
  This elementary algorithm leaves us with two conclusions. On the one hand, it is rather straightforward to calculate the expected value of the discrepancy function for a given box [0,x] × [0,y]. On the other hand, it is incredibly tedious to do so. While this calculation can be solved algorithmically in a straightforward fashion, there is little hope to compute the expectation analytically since we have a different set of success probabilities for each box generated by a vector (x,y).
  
  
  
  
  
  
  
  
  
  
   §.§ Numerical results
   
  In this final section, we present the results of three different sets of experiments. In the first two experiments we generate many instances of stratified point sets for a given fixed partition, calculate the _2-discrepancy of each point set and approximate the expected discrepancy of the partition by the mean of the experiment.
  In the final experiment, we calculate and compare the star discrepancy of different point sets.
  
  We use Warnock's formula <cit.> as presented in <cit.> to calculate the _2-discrepancy of a given point set. That is, for any point set ={ x_0, …, x_N-1}∈ [0,1]^d we have
  
      _2()^2 = 1/3^d - 2/N∑_n=0^N-1∏_i=0^d 1-x_n,i^2/2 + 1/N^2∑_m,n=0^N-1∏_i=0^d min(1-x_m,i, 1-x_n,i),
  
  in which x_n,i is the i-th component of the x_n. 
  We refer to <cit.> for quick implementations of this formula.
  
  First, we present a numerical observation which could, in principle, be proven along the same lines as Lemma <ref>. However, given the number of case distinctions such an analysis – based on our elementary method – would require, we only provide numerical evidence and state the result as a conjecture.
  
  
  There exists 𝐯=(v_1,v_2,v_3) with v_1< v_2< v_3 in [0,√(2)]  such that 
  
      _2^2(__𝐯^(4)) <  _2^2(_jit4) = 0.01909…
  
  
  
  We obtained various instances of partitions that seem to improve the classical jittered sampling by perturbing the three values of the vector (√(2)/4)(1,2,3). We obtained the best numerical results for the three points
  
      v_1^∗=√(2)/4 + 0.08,    v_2^∗=√(2)/2 + 0.11,    v_3^∗=3√(2)/4 - 0.02;
  
  see Fig. <ref>.
  We simulated 10^6 instances of stratified sets for this particular partition and calculated the discrepancy in each case with the formula of Warnock. Independently, we used our algorithm to estimate the expected discrepancy using 10^4 many grid points. 
  Both methods indicate that the first digits after the decimal points of the expected discrepancy are
  
      0.0188…,
  
  which would be clearly better than the mean discrepancy of jittered sampling. 
  
  
  Next, we use Warnock's formula to empirically study the discrepancy of different point sets and constructions; i.e. for given N we generate 500 samples and calculate the ℒ_2-discrepancy for each of these samples using Warnock's formula. The empirical mean of this sample approximates the expected value of the discrepancy.
  We collect our numerical results in Table <ref>.
  
  Our numerical results suggest that the expected discrepancy of partitions _∗^(N) is about a factor 2 smaller than the expected discrepancy of a set of random points:
   We conjecture that
  
      lim_N →∞_2(_N)^2/_2(__∗^(N))^2 = 2.
  
  
  
  
  
  
  
  
  In our final experiment, we use an implementation of the Dobkin-Eppstein-Mitchell algorithm <cit.>  for the computation of the star discrepancy which was provided by Magnus Wahlström; for details on the implementation we refer to <cit.>.
  This experiment relates to our comments in Section <ref> and shows that our partitions also seem to generate point sets that have a smaller expected star discrepancy than sets of N i.i.d. uniformly random points. We leave the generalisation of the partitions _∗^(N) for future research. It is in principle straightforward, but a bit technical and thus beyond the scope of this final proof-of-concept numerical experiment.
  
  
  
  
  
  
  20
  
  
  
  
  aisti C. Aistleitner, Covering numbers, dyadic chaining and discrepancy. J. Complexity 27 (2011), no. 6, 531-540. 
  
  chen1 W.W.L. Chen, M.M.Skriganov, Explicit constructions in the classical mean squares problem in irregularity of point distribution. J. Reine Angew. Math., 545 (2002), 67–95.
  
  cranley R. Cranley, T.N.L. Patterson, Randomization of number theoretic methods for multiple integration. SIAM Journal on Numerical Analysis 13(6) (1976), 904–914.
  
  volcic2 F. Chersi, A. Volčič, λ-equidistributed sequences of partitions and a theorem of the De Bruijn-Post type. Ann. Mat. Pura Appl., 162 (1992), 23–32.
  
  dave H. Davenport, Note on irregularities of distribution. Mathematika 3 (1956),131–135.
  
  dob D. P. Dobkin, D. Eppstein and D. P. Mitchell, Computing the Discrepancy with Applications to Supersampling Patterns, ACM Trans. Graph. (TOG) (15) (1996), no. 4, 354–376.
  
  persi P. Diaconis, The distribution of leading digits and uniform distribution mod 1, Ann. Probab. 5(1) (1977), 72–81.
  
  DP J. Dick, F. Pillichshammer, Digital Nets and Sequences, Cambridge Univ. Press, Cambridge, England, 2010.
  
  DP2 J. Dick, F. Pillichshammer, Explicit constructions of point sets and sequences with low discrepancy, Kritzer, Peter (ed.) et al., Uniform distribution and quasi-Monte Carlo methods. Discrepancy, integration and applications. Radon Series on Computational and Applied Mathematics 15, 63-86 (2014).
  
  DP3 J. Dick and F. Pillichshammer, Optimal ℒ_2-discrepancy bounds for higher order digital sequences over the finite field _2, Acta Arith. 162, No. 1 (2014), 65–99.
  
  dorr B. Doerr, A lower bound for the discrepancy of a random point set. J. Complexity 30 (2014), no. 1, 16–20. 
  
  doerr C. Doerr, M. Gnewuch and M. Wahlström, Calculation of Discrepancy Measures and Applications, Chen, William (ed.) et al., A panorama of discrepancy theory. Springer, Lecture Notes in Mathematics 2107, 621-678 (2014). 
  
  ecu2 P. L'Ecuyer, C. Lemieux, Variance reduction via lattice rules. Management Science 46(9) (2000), 1214–1235.
  
  ecu1 P. L'Ecuyer, Randomized Quasi-Monte Carlo: An Introduction for Practitioners. Owen, Art B. (ed.) et al., Monte Carlo and quasi-Monte Carlo methods, MCQMC 2016. Proceedings of the 12th international conference on Monte Carlo and quasi-Monte Carlo methods in scientific computing, Stanford, CA, August 14–19, 2016. Springer Proc. Math. Stat. 241, 29–52 (2018).
  
  Federer H. Federer, Curvature measures, Trans. Am. Math. Soc. 93 (1959), 418–491.
  
  heinrich1 K. Frank and S. Heinrich, Computing discrepancies of Smolyak quadrature rules, J. Complexity 12 (1996), 287–314.
  
  gnew M. Gnewuch,
  Bracketing numbers for axis-parallel boxes and applications to geometric discrepancy. 
  J. Complexity 24 (2008), no. 2, 154-172. 
  
  gnewuch M. Gnewuch, H. Pasing, C. Weiß, A Generalized Faulhaber Inequality, Improved Bracketing Covers, and Applications to Discrepancy, arXiv:2010.11479 
  
  haber S. Haber, A modified Monte Carlo quadrature. Math. Comput. 19 (1966), 361–368.
  
  
  
  heinrich2 S. Heinrich, Efficient algorithms for computing the L2 discrepancy. Math. Comput. 65 (1996), 1621–1633.
  
  hnww S. Heinrich, E. Novak, G. Wasilkowski and H. Wozniakowski, 
  The inverse of the star-discrepancy depends linearly on the dimension. 
  Acta Arith. 96 (2001), no. 3, 279–302. 
  
  hoeffding W. Hoeffding, On the distribution of the number of successes in independent trials. 
      Ann. Math. Statist. 27 (1956), 713–721.
  
  
  HuTaylor T-C Hu and R.L. Taylor, On the strong law for arrays and for the bootstrap mean and variance. Internat. J. Math. & Math. Sci., 20 (1997),  375–382. 
  
  kakutani S. Kakutani, A problem on equidistribution on the unit interval [0,1]. Measure Theory, Oberwolfach 1975, vol. 541, 369–375. Springer LNM (1976).
  
  klenke A. Klenke, Probability Theory – A Comprehensive Course. 2nd edn. Springer, Heidelberg (2013). 
  
  
  mat J. Matoušek, Geometric discrepancy. An illustrated guide. Algorithms and Combinatorics, vol 18. Springer, Berlin, Germany, 1999.
  
  nied H. Niederreiter, R. F. Tichy and G. Turnwald, An inequality for differences of distribution functions. Arch. Math. 54 (1990), 166–172.
  
  
  stefan1 F. Pausinger, S. Steinerberger, On the discrepancy of jittered sampling, J. Complexity 33 (2016), 199–216.
  stefan2 F. Pausinger, M. Rachh, S. Steinerberger, Optimal jittered sampling for two points in the unit square, Statist. Probab. Lett. 132 (2018), 55–61.
  
  owen1 A.B. Owen, Monte Carlo variance of scrambled equidistribution quadrature. SIAM Journal on Numerical Analysis 34(5) (1997), 1884–1910.
  owen2 A.B. Owen, Scrambled net variance for integrals of smooth functions. Annals of Statistics 25(4) (1997), 1541–1562.
  
  owen3 A.B. Owen, Variance with alternative scramblings of digital nets. ACM Transactions on Modeling and Computer Simulation 13(4) (2003), 363–378.
  
  
  
  pyke R. Pyke, W. R. van Zwet, Weak convergence results for the Kakutani interval splitting procedure, Ann. Probab. 32:1A (2004), 380–423.
  
  roth K. F. Roth, On irregularities of distribution. Mathematika 1 (1954), 73–79.
  
  schneider R. Schneider, Convex Bodies: The Brunn-Minkowsky Theory. 2nd edn. Cambridge Univ. Press, Cambridge, 2004. 
  
  SW R. Schneider, W. Weil, Stochastic and Integral Geometry. Springer, Berlin, 2008. 
  
  Thaele Ch. Thäle, 50 years sets with positive reacha survey. 
  Surv. Math.Appl. 3 (2008), 123–165.
  
  thompson 
   S.K. Thompson, Sampling. 3rd edn.  Wiley Series in Probability and Statistics
    Wiley, New York, 2012. 
  
  volcic A. Volčič, A generalization of Kakutani's splitting procedure. Ann. Mat. Pura Appl. 190:1 (2011), 45-54.
  
  
  
  warnock T. T. Warnock, Computational investigations of low discrepancy point sets. In S. K. Zaremba (editor), Applications of number theory to numerical analysis, pages 319–343. Academic Press, 1972.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § APPENDIX
  
  
  
  
   §.§ Appendix A – Uniform distribution and equidistribution of partitions
  
  
  The well-known fact that a sequence of Monte Carlo samples (_N) is almost surely uniformly distributed follows 
  directly from the strong law of large numbers. A corresponding statement for the sequences in Definition <ref> 
  is based on the strong law of large numbers for triangular arrays. 
  Note that it  does not require that the sampling points for different partitions ^(N) are independent. This is why we do not introduce this assumption, although it is typically satisfied in the applications we have in mind. 
  
  
  
  	Consider a sequence of partitions {^(N)}_N≥ 1 of [0,1]^d, with ^(N)=(Ω_1^(N),…, Ω_N^(N)) consisting of Lebesgue-sets with positive content and let ^(N)=(_1^(N),…, _N^(N)) be the vector of stratified sampling points based on ^(N). 
  	Then the triangular array =(^(N))_N∈ is almost surely uniformly distributed if and only if 
  	
      lim_N→∞∑_i=1^N 1/N|Ω_i^(N)||Ω_i^(N)∩ [ x,  y[|
      		=|[ x,  y[|
  
  	for all cubes [ x,  y[ ⊂ [0,1]^d.  
  
  
  	
  	For x,  y∈ [0,1]^d with each component of  at least as large as the corresponding component of , consider the  stochastic variables 
  	Y_i^(N)=1_[ x,  y[(^(N)_i), where ^(N)_1,…,_N^(N) is the stratified sample based on the partition ^(N). The Y's  are row-wise independent random variables with uniformly bounded variances and we have
  	
      a_N(,):=1/N∑_i=1^N Y_i^(N)=∑_i=1^N 1/N|Ω_i^(N)||Ω_i^(N)∩ [ x,  y[|.
  
  	Therefore, the strong law of large numbers for triangular arrays (see 
  	<cit.> with p=1, ψ(t)=t^2 and a_n=n) implies 
  	
      #(_^(N)∩ [ x,  y[)/N-a_N(,) → 0
  
  	almost surely as N→∞. 
  	
  	If assumption (<ref>) holds, we have lim_N→∞ a_N(,)=|[ x,  y[|, so (<ref>) implies  (<ref>) for almost every realization. Hence, the stratified sample points form almost surely a uniformly distributed triangular array. 
  	
  	If (<ref>)  is violated, there must be x,  y∈ [0,1]^d such that the limit in (<ref>) does not exist or is different from |[ x,  y[|. In any case, there is a subsequence (a_N'(,)) of (a_N(,)) such that a_N'(,)→ a |[ x,  y[| as N'→∞, and  (<ref>) shows that (<ref>) cannot hold for almost every realization. 
  	
  	Concluding,  the triangular array is uniformly distributed if and only if  (<ref>) holds 
  	for all rectangular sets [ x,  y[∈ [0,1]^d. 
  
  
  
  In particular,  if {^(N)}_N≥ 1 is a sequence of finite partitions of the unit cube such that all partitions are equivolume, then |Ω_i^(N)|=1/N, so (<ref>) is satisfied even without taking the limit. As a consequence, sequences of equivolume partitions are uniformly distributed; this is implication (d) in Fig. <ref>. 
  
  
  
  We conclude this section with a comparison of Definition <ref> with the notion of equidistributed partitions from  <cit.>, which we now recall using our notation.  
  A sequence {^(N)}_N≥ 1 of finite partitions of the unit cube is called equidistributed 
  if for any choice of  _i^(N)∈Ω_i^(N), i=1,…,N, the triangular array =(_1^(N),…,_N^(N))_N∈ is uniformly distributed (Definition <ref>). Actually, in <cit.> this notion is introduced and exploited in larger generality, replacing the unit cube with a general separable metric space. 
  
  
  In Fig. <ref>, we outline the connections between different notions for partitions, where a sequence {^(N)}_N≥ 1 of finite partitions of [0,1]^d is said to have vanishing average diameter if the average diameter 
  	
      δ_N=1/N∑_i=1^N diam (Ω_i^(N))
  
  	converges to zero as N→∞. We also use the set 
  	
      I_B(N)={i∈{1,…,N}: Ω_i^(N)⊂ B},
  
  	of all indices i for which Ω_i^(N) is completely contained in B⊂ [0,1]^2, and the condition 
  	
      lim_N→∞1/N#I_[ x,  y[(N)
      		=|[ x,  y[|,   for all [ x,  y[⊂ [0,1]^d,
  
  	stating that asymptotically the `correct' proportion of sets are completely contained in  [ x,  y[.  
  	That an equi­volume sequence of partitions with vanishing maximal diameter is equidistributed was shown in 
  	<cit.>, but our implication Fig. <ref>(a) is stronger, as it allows for `large' or `elongated' sets, as long as the proportion of such sets within the N sets of the partitions goes to zero as N→∞. 
  
  	
  		The implications in Fig. <ref> hold. Apart from the equivalence in (b) none of the implications can be reversed.	
  
  
  	To show (a) fix [ x,  y[ ⊂ [0,1]^d,
  		assume that the sets of {^(N)}_N≥ 1 are equivolume, and observe that 
  		in this case 
  		
      1/N#I_[ x,  y[(N)=| ⋃_i∈ I_
      				[ x,  y[(N)Ω^(N)_i|≤| ⋃_i=1^N(Ω^(N)_i∩ [ x,  y[)|=|[ x,  y[|.
  
  		On the other hand, if 
  		
      T_B(N)={i∈{1,…,N}: Ω_i^(N)∩ B∅,Ω_i^(N)∩ B^c∅},
  
  		describes the sets of the partition hitting both,  
  		B⊂ [0,1]^d and its complement, we have 
  		
      1/N#I_[ x,  y[(N)+1/N#T_[ x,  y[(N)
      		=| ⋃_i∈ I_[ x,  y[(N)∪ T_[ x,  y[(N)Ω^(N)_i|
      		≥ |[ x,  y[|,
  
  		 so it is enough to show that 
  		
      1/N#T_[ x,  y[(N)→ 0
  
  		as N→∞.  To show (<ref>) let ε>0 and α>0 be given. 
  		Markov's inequality and the assumption of vanishing average diameter yield the existence of  N_0∈ such that 
  		
      1/N#{i:diam (Ω_i^(N))>α}≤δ_N/α≤ε/2
  
  		for all N≥ N_0. Hence, using again the equivolume property,  
  		
      1/N#T_[ x,  y[(N)≤1/N#{i:diam (Ω_i^(N))>α}+
      		|( [ x,  y[)_α∩ [0,1]^d|≤ε/2+4dα,
  
  		where ( [ x,  y[)_α is the set of all points  with distance at most α from the boundary of [ x,  y[. Choosing α=ε/(8d) implies (<ref>).  
  	
  	
  	To show the equivalence in (b) we note that (<ref>) implies (<ref>). In fact, if  [ x,  y[ ⊂ [0,1]^d is fixed, the set B=[0,1[^d∖ [ x,  y[ can be written as disjoint union of at most k=3^d-1 half-open rectangles B_1,…,B_k, so 
  		
      1/N#T_[ x,  y[(N)≤1/N (N-#I_[ x,  y[(N)-∑_j=1^k #I_B_j(N))→ 0,
  
  		by applying   (<ref>) to all 3^d rectangles.
  	
  		If _N={_1,…,_N} satisfies _i∈Ω_i^(N) but is otherwise arbitrary,  then 
  		
      I_[,](N)≤#{i: _i∈ [0,]}≤ I_[,](N)+T_[,](N).
  
  		and there are sets for which either of the two inequalities becomes an equality. 
  		One can thus conclude that (<ref>) (together with its consequence, (<ref>)) implies equidistribution of {^(N)}_N≥ 1. Conversely, assuming equidistribution, we may choose _N such that there is equality on the left side of (<ref>) for every N, and  (<ref>)  follows.
  	
  	
  	Implication (c) is trivial as a stratified sample point lies a.s. in the partition set it is associated to.
  	Implication (d) has already been shown and (e) is trivial. 
  	None of the above implications can be reversed.That implication (a) cannot be reversed is clear, as one just takes a sequence satifying (<ref>) and changes two sets appropriately to assure that the sequence is not equivolume. Implication (d) can be treated in a similar way. That implication (c) cannot be reversed can be seen by means of vertical strip partitions (^(N)_vert) in
  	Fig.  <ref> (left). For all N these partitions are equivolume and thus the stratified sample based on them is a.s. uniformly distributed. However, this sequence of partitions is not equidistributed, as the limit in (<ref>) is always zero. This example also shows that implication (e) cannot be reversed.
  
  
  
  
  
   §.§ Appendix B – Proof of Lemma <ref>
  
  In this part of the appendix we present the details we omitted in the proof of Lemma <ref> in Section <ref>. Lemma <ref> can then be shown along the very same lines with slightly different probabilities and boundaries for the six subcases; see Fig. <ref>.
  
  
  
  We set A=v_1 √(2), B=v_2 √(2)-1 such that for given 1/2 ≤ A ≤ 1 we have 2A-1 ≤ B ≤ A.
  
  Recall that we 
  	
  	
  	need to evaluate (<ref>) with the notation introduced just before this formula was stated.
  The probabilities in each of the  cases below are obtained from elementary calculations.
  
  
  
  
    
  Case I The case (x,y) ∈Ω_1=S_1 was already discussed in the main text.
  
  
  
    
  Case II 
  Let (x,y) ∈ S_2; i.e. 0≤ x ≤ B and A ≤ y ≤ 1 as well as B≤ x ≤ A and A ≤ y ≤ 1+B-x. Then q_3(x,y)=0,
  
      q_1(x,y)=2 (A x-x^2/2)/A^2    and    
      q_2(x,y)=x (y-A)+x^2/2/-A^2/2-1/2 (1-B)^2+1.
   
  Hence, we get
  
      f_2   (x,y)= 1/9 A^2 (A^2+(B-2)
         B-1)×
         ( x^3 (-2 y (-6 A^2-3 (B-2) B+1)-8A) 
         +x^2 (3 A^2 y^2 (3 A^2+3 (B-2) B+1)-4 A y (6 A^2+3 (B-2) B-1)+6 A^2+(2-B) B+1) 
         +x (4 A^3-2 A^2 y+2 A (B-2) B-2 A) +2 x^4b).
  
  
  
  
  
  
  
    
  Case III 
  Let (x,y) ∈ S_3; i.e. 0 ≤ x ≤ A/2 and A-x ≤ y ≤ A as well as A/2 ≤ x ≤ A and x ≤ y ≤ A. Then q_3(x,y)=0 and 
  
      q_1(x,y)=2 (x y-1/2 (-A+x+y)^2)/A^2    and    
      q_2(x,y)=-(-A+x+y)^2/A^2+(B-2) B-1.
   
  Hence, we get
  
      f_3   (x,y)= 1/9 A^2 (A^2+(B-2) B-1)×
         ( 2 x^4 + x^3 (12 A^2 y-8 A+6 (B-2) B y-2
         y) . 
         +x (12 A^4 y-24
         A^3 y^2-4 A^3+6 A^2 B^2 y-12 A^2 B y+12 A^2 y^3+12A^2 y-12 A (B-2) B y^2  . 
             . +2 A (B-2) B-4 A y^2-2 A+6 (B-2) B y^3-2 y^3) 
          +x^2 (9 A^4 y^2-24 A^3 y+3 A^2 (3 (B-2) B+1) y^2+10 A^2-12 A (B-2) B y . 
             . -4 A y+(2-B) B+4 y^2+1) 
           +2 A (B-2) B y-8 A y^3-2 A  y-(B-2) B y^2+2 y^4+y^2 +2 A^2 B 
             . +10 A^2 y^2+A^2 -4 A^3 y-A^2 B^2 ).
  
  
  
  
  
  
  
    
  Case IV 
  Let (x,y) ∈ S_4; i.e. A ≤ x ≤ 1+B-A and A ≤ y ≤ 1+B-x. Then q_1(x,y)=1, q_3(x,y)=0 and
  
      q_2(x,y)= x y-x^2/2/-A^2/2-1/2 (1-B)^2+1.
   
  Hence, we get
  
      f_4(x,y)= (-1 - 2 B (1 - 3 x y)^2 + B^2 (1 - 3 x y)^2 +  x^2 (4 + 3 y (-4 x + y + 3 A^2 y)))/9 (-1 + A^2 + (-2 + B) B).
  
  
  
  
  
  
  
    
  Case V 
  Let (x,y) ∈ S_5; i.e. B ≤ x ≤ A and 1+B-x ≤ y ≤ 1. Then q_2(x,y)=q_3(x,y)=0 and q_1(x,y)=2xy/A^2. Hence, we get
  
      f_5   (x,y)= 1/9 A^2 (-1 + B)^2 (-1 + A^2 + (-2 + B) B)×
         ( A^4 B^2-2 A^4 B y+2 A^4 B+A^4 y^2-2 A^4 y+A^4+2 A^2 (B+1)^2 (2 B^2+1)
       . 
              -8 A^2 B y^3+2 A^2 (B (7 B+10)+6) y^2-4 A^2 B (B (3 B+5)+4) y+2 A^2 y^4-8 A^2
         y^3-8 A^2 y 
          +x^4 (-2 A^2-2 B^2+4 B+2) 
            +x^3 (-6 A^4 y+8 A^3+8 A^2 y-8 A+6
         B^4 y+B^3 (8-24 y)-2 B^2 (4-10 y)+8 B (y-2)-2
         y) 
            +x^2 (9 A^4 (B-2) B y^2+12 A^4 B y-3
         A^4 y^2+12 A^4 y+A^4-24 A^3 (B-2) B y-16 A^3 B . 
               -8 A^3 y-16 A^3+2 A^2 (8 B^2+7)+3 A^2 ((B-2)
         B-1) (3 (B-2) B-1) y^2 
               +8 A^2 (B+1) (3 (B-2) B-1)y-12 A B^4 y+8 A B^3 (6 y-2)+8 A B^2 (2-5 y) 
             . -16 A B
         (y-2)+4 A y-5 B^4+8 B^3 y+4 B^3-4 B^2 y (y+2)+8
         B^2+8 B (y-1)^2+1) 
            +x (-6 A^4 B^2 y+12 A^4
         B y^2-12 A^4 B y-2 A^4 B-6 A^4 y^3+12 A^4 y^2-4 A^4
         y-2 A^4+12 A^3 B^2  . 
               -16 A^3 B y+8 A^3 B+8 A^3 y^2-16
         A^3 y+12 A^3-2 A^2 B (6 B^3-29 B-30) y 
              -12 A^2 (B-2) B y^3 +8 A^2 (B+1) (3 (B-2) B-2) y^2-4 A^2
         (B+1) (B (3 B+2)+2) 
              +4 A^2 y^3+18 A^2 y+10 A B^4-16 A
         B^3 y-8 A B^3+8 A B^2 y (y+2)
         
             . .-16 A B^2-16 A B
         (y-1)^2-2 A) ).
  
  
  
  
  
  
  
  
  
    
  Case VI 
  Let (x,y) ∈ S_6; i.e. A ≤ x ≤ (B+1)/2 and 1+B-x ≤ y ≤ 1 as well as (B+1)/2 ≤ x ≤ 1 and x ≤ y ≤ 1. Then q_1(x,y)=1 and 
  
      q_2(x,y)=-A^2/2-1/2 (-B+x+y-1)^2+x
         y/-A^2/2-1/2 (1-B)^2+1    and     q_3(x,y)=(-B+x+y-1)^2/(1-B)^2.
   
  Hence, we get
  
      f_6    (x,y)= 1/9 (B-1)^2 (A^2+(B-2) B-1)×
         ( A^4 B^2-2 A^4 B y+2 A^4 B+A^4 y^2-2 A^4 y+A^4 +2 A^2 (B+1)^2 (2 B^2+1) . 
            -8 A^2 B y^3+2 A^2 (B (7 B+10)+6) y^2-4 A^2 B (B (3 B+5)+4) y+2 A^2 y^4-8 A^2 y^3-8 A^2 y 
         +x^4
         (-2 A^2-2 B^2+4 B+2) 
            +x^3 (-6 A^4 y+8 A^3+8 A^2 y-8 A+6
         B^4 y+B^3 (8-24 y)-2 B^2 (4-10 y)+8 B (y-2)-2
         y) 
          +x^2 (9 A^4 (B-2) B y^2+12 A^4 B y-3
         A^4 y^2+12 A^4 y+A^4-24 A^3 (B-2) B y-16 A^3 B . 
               -8 A^3y-16 A^3+2 A^2 (8 B^2+7)+3 A^2 ((B-2)
         B-1) (3 (B-2) B-1) y^2
              +8 A^2 (B+1) (3 (B-2) B-1)y -12 A B^4 y+8 A B^3 (6 y-2)+8 A B^2 (2-5 y) 
             . -16 A B
         (y-2)+4 A y-5 B^4+8 B^3 y+4 B^3-4 B^2 y (y+2)+8
         B^2+8 B (y-1)^2+1) 
          +x (-6 A^4 B^2 y+12 A^4
         B y^2-12 A^4 B y-2 A^4 B-6 A^4 y^3+12 A^4 y^2-4 A^4
         y-2 A^4+12 A^3 B^2 . 
              -16 A^3 B y+8 A^3 B+8 A^3 y^2-16
         A^3 y+12 A^3-2 A^2 B (6 B^3-29 B-30) y 
              -12  A^2 (B-2) B y^3+8 A^2 (B+1) (3 (B-2) B-2) y^2-4 A^2
         (B+1) (B (3 B+2)+2) 
               +4 A^2 y^3+18 A^2 y+10 A B^4-16 A
         B^3 y-8 A B^3+8 A B^2 y (y+2)
         
             . .
         -16 A B^2-16 A B
         (y-1)^2-2 A) ).
  
  
  
  
  
  
  Inserting the explicit expressions of the functions f_1, …, f_6 into (<ref>) 
  gives us the following rational function in A and B which we can minimize over 1/2≤ A ≤ 1 and 2A-1 ≤ B ≤ A in order to obtain the two parameters A and B that generate the partition with the smallest expected discrepancy in this family: 
  
      _2^2(__v_1,v_2^(3) ) =
      1/25920 A^2 (A^2+(B-2) B-1)×
         . (480 A^8+6912 A^7 (B+1)+16 A^6 (-527 + 330 B + 195 B^2) . 
         . -384 A^5 (13 + 67 B + 87 B^2 + 33 B^3)+ . 
         . 12 A^4 (829 + 1952 B + 2720 B^2 + 1228 B^3 + 431 B^4) . 
         . +48 A^3 (-54 - 181 B - 204 B^2 + 23 B^3 + 120 B^4 + 24 B^5) . 
         . -4 A^2 (118 + 104 B + 345 B^2 + 776 B^3 - 728 B^4 + 1368 B^5 + 393 B^6). 
         . +48 A (-10 + 7 B + 76 B^2 - 3 B^3 - 140 B^4 - 42 B^5 + 36 B^6 + 12 B^7). 
         . -(56 B + 774 B^2 + 556 B^3 - 1364 B^4 - 1728 B^5 - 70 B^6 + 396 B^7 + 
       99 B^8)).
  
  This function can be minimised using a standard computer algebra system.
  The minimum of this function is 0.0268044 for the parameter values A= 0.7512174 and B= 0.513013; see Fig. <ref>. 
  These parameter values correspond to v_1 = 0.5311… and v_2 = 1.0698….
  
  
  
  
  