A rejection technique for sampling from log-concave multivariate distributions.
Different universal methods (also called automatic or black-box methods) have been suggested for sampling form univariate log-concave distributions. The descriptioon of a suitable universal generator for multivariate distributions in arbitrary dimensions has not been published up to now. The new algorithm is based on the method of transformed density rejection. To construct a hat function for the rejection algorithm the multivariate density is transformed by a proper transformation T into a concave function (in the case of log-concave density log(x).)  Then it is possible to construct   a dominating function by taking the minimum of serveral tangent hyperplanes that are transformed back by  T-1 into the  original scale. The domains of different pieces of the hat function  are polyhedra in the multivariate case. Although this method can be shown to work, it is too slow and complicated in higher dimensions. In this article we split the  Rn into simple cones. The hat function is constructed piecewise on each of the cones by tangent hyperplanes. The resulting function is no longer continuous and the rejection constant is bounded from below but the setup and the generation remains quite fast in higher dimensions; for example, 8. The article describes the details of how this main  idea can be used to construct algorithm TDRMV that generates random tuples from a multivariate log-concave distribution  with a computable density. Although the developed  algorithm is  not a real black box method it is adjustable for a large class of log-concave densities.
Introduction
For the univariate case there is a large literature on generation methods for
standard distributions (see e.g. [Dev86] and [Dag88]) and in the last years some
papers appeared on universal (or black-box) methods (see [Dev86, chapter VII],
[GW92], [Ahr95], [H-or95a], [HD94] and [ES97]); these are algorithms that can
generate random variates from a large family as long as some information (typ-
ically the mode and the density of the specific distribution) are available.
For the generation of variates from bivariate and multivariate distributions
papers are rare. Well known and discussed are only the generation of the multi-normal
and of the Wishart distribution (see e.g. [Dev86] and [Dag88]). Several
approaches to the problem of generating multivariate random tuples exist, but
these have some disadvantages:
ffl The multivariate extension of the ratio of uniforms methods as in [SV87]
or [WGS91]. This method can be reformulated as rejection from a small
family of table-mountain shaped multivariate distributions. This point of
view is not included in these two papers but it is useful as it clarifies the
question why the acceptance probability becomes poor for high correla-
tion. This disadvantage of the method is already mentioned in [WGS91].
The practical problem how to obtain the necessary multivariate rectangle
enclosing the region of acceptance for the ratio of uniforms method is not
discussed in [SV87] nor in [WGS91] and seems to be difficult for most
distributions.
ffl The conditional distribution method. It requires the knowledge of and the
ability to sample from the marginal and the conditional distributions (see
[Dev86, chapter XI.1.2]).
ffl The decomposition and rejection method. A majorizing function (also
called suggested for the multivariate rejection method is the
product of the marginal densities (in [Dag88]). It is not clear at all how
to obtain the necessary rejection constant ff.
ffl Development of new classes of multivariate distributions, which are easy
to generate. It is only necessary (and possible) to specify the marginal
distribution and the degree of dependence measured by some correlation
coefficient (see the monograph [Joh87]). This idea seems to be attractive
for most simulation practitioners interested in multivariate distributions
but it is no help if variates from a distribution with given density should
be generated.
Recently Devroye [Dev97] has developed algorithms for ortho-unimodal
densities. But this paper leaves the generation of log-convave distributions
as open problem.
ffl Sweep-plane methods for log-concave (and T-concave) distributions are
described recently in [H-or95b] for bivariate case and in [LH98] for the
multivariate case. These algorithms use the idea of transformed density
rejection which is presented in a first form in [Dev86, chapter VII.2.4] and
with a different set-up in [GW92].
To our knowledge these two algorithms are the only universal algorithms
in the literature for multivariate distributions with given densities. (In
[Dev86, chapter XI.1.3] it is even stressed that no general inequalities
for multivariate densities are available, a fact which makes the design of
black-box algorithms, that are similar to those developed in [Dev86] for
the univariate case, impossible.)
Although the algorithm in [LH98] works, it is very slow, since the domain of
the density f is decomposed in polyhedra. This is due to the construction of the
hat function, where we take the pointwise minimum of tangent hyperplanes. In
this paper we again use transformed density rejection and the sweep-plane technique
to derive a much more efficient algorithm. The main idea is to decompose
the domain of the density in cones first and then compute tangent hyperplanes
in this cones. The resulting hat function is not continuous any more and the
rejection constant is bounded from below, but the setup as well as the sampling
from the hat function is much faster than in the original algorithm.
Section 2 explains the method and gives all necessary mathematical formu-
lae. Section 3 provides all details of the algorithm. Section 4 discusses how to
improve and extend the main idea of the algorithm (e.g. to T-concave distribu-
tions, bounded domain) and section 5 reports the computational experience we
have had with the new algorithm.
2 The method
2.1 Transformed density rejection
Density. We are given a multivariate distribution with differentiable density
function
To simplify the development of our method we assume
In x4 we extend the algorithm so that these requirements can be
dropped.
Transformation. To design an universal algorithm utilizing the rejection
method it is necessary to find an automatic way to construct a hat function
for a given density. Transformed density rejection introduced under a different
name in [GW92] and generalized in [H-or95a] is based on the idea that the density
f is transformed by a monotone T (e.g. T in such a way that
(see [H-or95a]):
concave (we then say "f is T -concave");
differentiable and T 0 (x) ? 0, which implies T \Gamma1 exists; and
(T4) the volume under the hat is finite.
Hat. It is then easy to construct a hat ~ h(x) for ~
f(x) as the minimum of N
tangents. Since ~
f(x) is concave we clearly have ~
Transforming ~ h(x) back into the original scale we get
majorizing function or hat for f , i.e. with f(x) - h(x). Figure 1 illustrates
the situation for the univariate case by means of the normal distribution and
the transformation T log(x). The left hand side shows the transformed
Figure

1: hat function for univariate normal density
density with three tangents. The right hand side shows the density function
with the resulting hat. (The dashed lines are simple lower bounds for the density
called squeezes in random variate generation. Their use reduces the number of
evaluations of f . Especially if the number of touching points is large and the
evaluation of f is slow the acceleration gained by the squeezes can be enormous.)
Rejection. The basic form of the multivariate rejection method is given by
algorithm Rejection().
Algorithm 1 Rejection()
1: Set-up: Construct a hat-function h(x).
2: Generate a random tuple proportional to
h(X) and a uniform random number U .
3: If Uh(X) - f(X) return X else go to 2.
The main idea of this paper is to extend transformed density rejection as
described in [H-or95a] to the multivariate case.
2.2 Construction of a hat function
Tangents. Let p i be points in D ' R n . In the multivariate case the tangents
of the transformed density ~
f(x) at p i are the hyperplanes given by
where h\Delta; \Deltai denotes the scalar product.
Polyhedra. In [LH98] a hat function h(x) is constructed by the pointwise
minimum of these tangents. We have
The domains in which a particular tangent ' i (x) determines the hat function
are simple convex polyhedra P i , which may be bounded or not (for details
about convex polyhedra see [Gr-u67, Zie95]). Then a sweep-plane technique for
generating random tuples in such a polyhedron with density proportional to
To avoid lots of indices we write p, '(x) and P without the index i if there
is no risk of confusion.
A sweep-plane algorithm. Let
r ~
kr ~
if r ~
choose any g with denotes the 2-norm.)
For a given x let We denote the hyperplane perpendicular to g
through x by
fy
and its intersection with the polytope P with
depends on x only; thus we write F (x), if there is no risk of confusion.) Q(x)
again is a convex simple polyhedra. Now we can move this sweep-plane F (x)
through the domain P by varying x. Figure 2 illustrates the situation.
As can easily be seen from (2), (4) and (5), T \Gamma1 ('(x)) is constant on Q(x)
for every x. Let
Then the hat function in P is given by
where again We find for the marginal density function of the hat
Z

Figure

2: sweep-plane F (x)
where integration is done over F (x). A(x) denotes the (n \Gamma 1)-dimensional
volume of Q(x). It exists if and only if Q(x) is bounded.
To compute A(x) let v j denote the vertices of P and v
assume that the polyhedron P is simple. Then let t v j
n be the n nonzero
vectors in the directions of the edges of P originated from v j , i.e. for each k and
every
by modifying the method in [Law91] we find
a
The coefficients are given by
Y
and
a
Notice that b (x)
equations (9) and (10)
does not hold if P is not simple. For details see [LH98].
The generation from h g
is not easy in general. But for log-concave or
-concave (see x4.8) densities f(x), h g
again is log-concave ([Pr'e73]) and T c -
concave ([LH98]), respectively.
Generate random tuples. For sampling from the "hat distribution" we first
need the volume below the hat in all the polyhedra P i and in the domain D.
We then choose one of these polytopes randomly with density proportional to
their volumes. By means of a proper univariate random number we sample from
marginal distribution hj g
and get a intersection Q(x) of P . At last we have to
sample from a uniform distribution on Q(x).
It can be shown (see [LH98]) that the algorithm works if
(1) the polyhedra P i are simple (see above),
(2) there exists a unique maximum of ' i (x) in P i (then ff \Gamma fi x is decreasing
and thus the volume below the hat is finite in unbounded polyhedra), and
is non-constant on every edge of P i (otherwise hg; t v j
and an edge t i and thus a
Adaptive rejection sampling. It is very hard to find optimal points for
constructing these tangents ' i (x). Thus these points must be chosen by adaptive
rejection sampling (see [GW92]). Adapted to our situation it works in the
following way: We start with the vertices of a regular simplex and add a
new construction point whenever a point is rejected until the maximum number
N of tangents is reached. The points of contact are thus chosen by a stochastic
algorithm and it is clear that the multivariate density of the distribution of
the next point for a new tangent is proportional to h(x) \Gamma f(x). Hence with
tending towards infinity the acceptance probability for a hat constructed in
such a way converges to 1 with probability 1. It is not difficult to show that the
expected volume below the hat is 1 +O(N \Gamma2=n ).
Problems. Using this method we run into several problems.
We have to compute the polyhedra every time we add a point.
What must be done, if the marginal distribution (8) does not exist in the
initial (usually not bounded) polyhedra P i , or if the volume below the hat
is infinite (Q i (x) not bounded, ff \Gamma fi x not decreasing)?
Moreover the polyhedra P i typically have many vertices. Therefore the algorithm
is slow and hard to implement because of the following effects.
\Gamma The computation of the polyhedra (setup) is very expensive.
\Gamma The marginal density (8) is expensive to compute. Since it is different
for every polyhedron P i (and for every density function f ), we have to
use a slow black box method (e.g. [GW92, H-or95a]) for sampling from the
marginal distribution even in the case of log-concave densities.
\Gamma Q(x) is not a simplex. Thus we have to use the (slow) recursive sweep-
plane algorithm as described in [LH98] for sampling from the uniform
distribution over a (simple) polytope.
2.3 Simple cones
A better idea is to choose the polyhedra first as simple as possible, i.e. we choose
cones. (We describe in x2.4 how to get such cones.)
A simple cone C (with its vertex in the origin) is an unbounded subset
spanned by n linearly independent vectors:
In opposition to the procedure described above we now have to choose a proper
point p in this cone C for constructing a tangent. In the whole cone the hat h
is then given by this tangent. The method itself remains the same.
Obviously the hat function is not continuous any more (because we first define
a decomposition of the domain and then compute the hat function over the
different parts. It cannot be made continuous by taking the pointwise minimum
of the tangents, since otherwise we cannot compute the marginal density h g
by
equation (8)). Moreover we have to choose one touching point in each part.
These disadvantages are negligible compared to the enormous speedup of the
setup and of the generation of random tuples with respect to this hat function.
Marginal density. The intersection Q(x) of the sweep plane F (x) with the
cone C is bounded if and only if F (x) cuts each of the sets f-t
x ? 0, i.e. if and only if hg; -t hence if and only if
We find for the volume A(x) in (9) of the intersection Q(x)
ae a x
where (again)
a
Y
Notice that A(x) does not exist if condition (13) is violated, whereas the right
hand side in (14) is defined
If the marginal density exists, i.e. (13) holds, then by (8) and (6) it is given
by
Volume. The volume below the hat function in a cone C is given by
Z 1h g
Z 1a x
Notice that g and thus a, ff and fi depend on the choice of p. Choosing an
arbitrary p may result in a very large volume below the hat and thus in a very
poor rejection constant.
Intersection of sweep-plane. Notice that the intersection Q(x) is always
a (n \Gamma 1)-simplex if condition (13) holds. Thus we can use the algorithm in
[Dev86] for sampling from uniform distribution on Q(x). The vertices
of Q(x) in R n are given by
uniformly [0; 1] random variates and U
We sort these variates such that U 0 - Un . Then we get a random point
in Q(x) by (see [Dev86, theorems XI.2.5 and V.2.1])
The choice of p. One of the main difficulties of the new approach is the choice
of the touching point p. In opposite to the first approach where the polyhedron
is build around the touching point, we now have to find such a point so that
holds. Moreover the volume below the hat function over the cone should
be as small as possible.
Searching for such a touching point in the whole cone C or in domain D
(the touching point needs not to be in C) with techniques for multidimensional
minimization is not very practicable. Firstly the evaluation of the the volume
HC in (17) for a given point p is expensive and its gradient with respect to p is
not given. Secondly the domain of HC is given by the set of points where (13)
holds.
Instead we suggest to choose a point in the center of C for a proper touching
point for our hat. Let -
be the barycenter of the spanning vectors.
Let a(s), ff(s) and fi(s) denote the corresponding parameters in (16) for t.
Then we choose by minimizing the function
Z
The domain DA of this function is given by all points, where kr ~
where A(x) exists, i.e. where condition (13). It can easily be
seen, that DA is an open subset of (0; 1).
To minimize j we can use standard methods, e.g. Brent's algorithm (see
e.g. [FMM77]). The main problem is to find DA . Although ~
f(x) is concave by
assumption, it is possible for a particular cone C that DA is a strict subset of
(0; 1) or even the empty set. Moreover it might not be connected. In general
only the following holds: Let (a; b) be a component of DA 6= ;. If If f 2 C 1 , i.e.
the gradient of f is continuous, then
lim
s&a
Roughly spoken, j is a U-shaped function on (a; b).
An essential part of the minimization is initial bracketing of the minimum,
i.e. finding three points s (a; b), such that j(s 1
This is necessary since the function term of j in (20) is also
defined for some s 62 DA (e.g. s ! 0). Using Brent's algorithm without initial
bracketing may (and occasionally does) result in e.g. a negative s.
Bracketing can be done by (1) search for a s 1 2 DA , and (2) use property
(21) and move towards a and b, respectively, to find an s 0 and an s 2 . (It is
obvious that we only find a local minimum of j by this procedure. But in all
the distributions we have tested, there is just one local minimum which therefore
is the global one.)
For the special case where hg(s); - ti does not depend on s (e.g. for all multivariate
normal distributions) DA either is (0; 1) or the empty set. It is then
possible to make similar considerations like that in [H-or95a, theorem 2.1] for
the one dimensional case. Adapted to the multivariate case it would state, that
for the optimal touching point p, f(p) is the same for every cone C.
Condition violated. Notice that DA even may be the empty set, i.e. condition
fails for all s 2 (0; 1). By the concavity of ~
f(x) we know, that
construction point p. Furthermore hg; pi is bounded from
below on every compact subset of the domain D of the density f . Therefore
there always exists a partition into simple cones with proper touching points
which satisfy (13), i.e. the domains DA are not empty for all cones C.
We even can have
2.4 Triangulation
For this new approach we need a partition of the R n into simple cones. We get
such a partition by triangulation of the unit sphere S n\Gamma1 . Each cone C is then
generated by a simplex \Delta ae S n\Gamma1 (triangle in S 2 , tetrahedron in S 3 , and so
\Deltag (22)
These simplices are uniquely determined by the vectors t
their vertices. (They are the the convex hull of these vertices in S n\Gamma1 .) It does
not matter that these cones are closed sets. The intersection of such cones might
not be empty but has measure zero.
For computing a in (15) we need the volumes of these simplices. To avoid DA
being the empty set, some of the cones have to be skinny. Furthermore to get
a good hat function, these simplices should have the same volume (if possible)
and they should be "regular", i.e. the distances from the center to the vertices
should be equal (or similar). Thus the triangulation should have the following
properties:
Recursive construction.
are easy computable for all simplices.
(C3) Edges of a simplex have equal length.
Although it is not possible to get such a triangulation for n - 3 we suggest an
algorithm which fulfils (C1) and (C2) and which "nearly" satisfies (C3).
Initial cones. We get the initial simplices as the convex hull in S n\Gamma1 of the
vectors
en (23)
where e i denotes the i-th unit vector in R n (i.e. a vector where the i-th component
is 1 and all others are 1g. As can easily be seen the
resulting partition of the R n is that of the arrangement of the hyperplanes
Hence we have 2 n initial cones.
Barycentric subdivision of edges. To get smaller cones we have to triangulate
these simplices. Standard triangulations of simplices which are used
for example in fixed-point computation (see e.g. [Tod76, Tod78]) are not appropriate
for our purpose. The number of simplices increases too fast for each
triangulation step. (In opposition to fixed point calculations, we have to keep
all simplices with all their parameters in the memory of the computer.)
Instead we use a barycentric subdivision of edges: Let t be the
vertices of a simplex \Delta. Then use the following algorithm.
(1) Find the longest edge
(2) Let
i.e. the barycenter of the edge projected to the sphere.
(3) Get two smaller simplices: Replace vertex t i by t new for the first simplex
and vertex t j by t new for the second one. We have
After making k of such triangulation steps in all initial cones we have 2 n+k
simplices.
This triangulation is more flexible. Whenever we have a cone C, where D a is
empty (or the algorithm does not find an s 2 D a ) we can split C and try again
to find a proper touching point in both new cones. This can be done until we
have found proper construction points for all cones of the partition (see end of
x2.3). In practice this procedure stops, if too many cones are necessary. (The
computer runs out of memory.)
Notice that it is not a good idea to use barycentric subdivision of the whole
simplex (instead of dividing the longest edge). This triangulation exhibits the
inefficient behavior of creating long, skinny simplices (see remark in [Tod76]).
"Oldest" edge. Finding the longest of the
edges of a simplex is very
expensive. An alternative approach is to use the "oldest" edge of a simplex.
The idea is the following:
(1) Enumerate the 2n vertices of the initial cones.
(2) Whenever a new vertex is created by barycentric subdivision, it gets the
next number.
(3) Edges are indexed by the tuple (i; j) of the number of the incident vertices,
such that i ! j.
We choose the edge with the lowest index with respect to the lexicographic
order (the "oldest" edge). This is just the pair of lowest indices of the
vertices of the simplex.
As can easily be seen, the "oldest" edge is (one of) the longest edge(s) for the
first steps. Unluckily this does not hold for all simplices in
following triangulation steps. (But it is at least not the shortest one.)
Computational experiences with several normal distributions for some dimensions
have show, that this idea speeds up the triangulation enormously
but has very little effect on the rejection constant.
Setup. The basic version of the setup algorithm is as follows:
1. Create initial cones.
2. Triangulate.
3. Find touching points p if possible (and necessary).
4. Triangulate every cone without proper touching point.
5. Goto 3 if cones without proper touching points exist, otherwise stop.
2.5 Problems
Although this procedure works for our tested distributions, an adaptation might
be necessary for a particular density function f .
(1) The searching algorithm for a proper touching point in x2.3 can be im-
proved. E.g. DA is either [0; 1) or the empty set if f is a normal distribution.
(2) There is no criterion how many triangulation steps are necessary or usefull
for an optimal rejection constant. Thus some tests with different numbers of
trianglation steps should be made with density f (see also x5).
(3) It is possible to triangulate each cone with a "bad" touching point. But
besides the case where no proper touching point can be found, some touching
points may lead to an enourmous volume below the hat function. So this case
should also be excluded and the corresponding cones should be triangulated.
A simple solution to this problem is that an upper bound Hmax for the
volumes HC is provided. Each cone with HC ? Hmax has to be triangulated
further. Such a bound can be found by some empirical tests with the given
density f .
Another way is to triangulate all initial cones first and then let Hmax be a
multiple (e.g. 10) of the 90th percentile of the HC of all created cones.
Problems might occur when the mode is on the boundary of the support
(Then we set
can be seen as a concave
An example for such a situation is when f(x) is a normal density on
a ball B and vanishes outside of B.
In such a case there exists a cone C such that f- t: - ? 0g does not intersect
suppf and the algorithm is in troubles. If C " we simply can
remove this cone. Otherwise an expensive search for a proper touching point is
necessary.
Restrictions. The above observations - besides the fact that no automatic
adaption is possible - are a drawback of the algorithm for its usage as black-box
algorithm. Nevertheless the algorithm is suitable for a large class of log-concave
densities and it is possible to include parameters into the code to adjust the
algorithm for a given density easily. Of course some tests might be necessary.
Besides, the algorithm does not produce wrong random points but simply does
not work, if no "good" touching points can be found for some cones C.
2.6 Log-concave densities
The transformation T
is concave, we say f is log-concave. We have T thus we find
for the marginal density function in (16) those of a gamma distribution with
shape parameters n and fi.
The volume below the hat for log-concave densities in a cone C is now given by
Z 1a x
To minimize this function it is best to use its logarithm:
For the normal distribution with density proportional to
we have ~
is the center of
the cone C with
1. Thus we simply find by (6)
Since a(s) does not depend on s we find for constant.
But even for the normal distribution with an arbitrary covariance matrix, this
function becomes much more complicated.
3 The algorithm
The algorithm tdrmv() consists of two main parts: the construction of a hat
function h(x) and the generation of random tuples X with density proportional
to this hat function. The first one is done by the subroutine setup(), the second
one by the routine sample().
Algorithm random tuple for given log-concave density
Input: density f
1: call Construct a hat-function h(x)
2: repeat
3: X / call sample(). =  Generate a random tuple X with density prop. to h(X).
4: Generate a uniform random number U .
5: until U \Delta h(X) - f(X).
return X.
To store h(x), we need a list of all cones C. For each of these cones we need
several data which we store in the object cone. Notice that the variables p, g,
ff, fi, a and HC depend on the choice the touching point p and thus on s. Some
of the parameters are only necessary for the setup.
object 1 cone
parameter variable definition
spanning vectors t
center of cone -
construction point
location of p s (setup)
sweep plane g see (4)
marginal density ff; fi see (6)
coefficient a see (15)
determinant of vectors
volume under hat HC ; H cum
C see (17) and (28)
Remark. To make the description of the algorithm more readable, some standard
techniques are not given in details.
3.1
The routine setup() consists of three parts: (H1) setup initial cones, (H2)
triangulation of the initial cones and (H3) calculation of parameters.
(H1) is simple (see x2.4). (H2) is done by subroutine split(). The main
problem in (H3) is how to find the parameter s (i.e. a proper construction
point). This is done by subroutine find(). Minimizing (29) is very expensive.
Notice that for a given s we have to compute all parameters that depend on s
before evaluating this function. Since it is not suitable to use the derivative of
this function, a good choice for finding the minimum is to use Brent's algorithm
(e.g. [FMM77]). To reduce the cost for finding a proper s, we do not minimize
for every cone. Instead we use the following procedure:
(1) Make some triangulation steps as described in x2.4.
(2) Compute s for every cone C.
(3) Continue with triangulation. When a cone is split by barycentric subdivision
of the corresponding simplex, both new cones inherit s from the old simplex.
Our computational experiences with various normal distributions show, that
the costs for setup reduces enormously without raising the rejection constant
too much. Using this procedure it might happen that s does not give a proper
touching point (or HC is too big; see end of x2.4) after finishing all
triangulation steps. Thus we have to check s for every cone and continue with
triangulation in some cones if necessary.
3.2 Sampling
The subroutine sample() consists of four parts: (S1) select a cone C, (S2) find
a random variate proportional to the marginal density h g
(27), (S3) generate a
uniform random tuple U on the standard simplex (i.e.
and compute tuple on the intersection Q(x) of the
sweeping plane with cone C. (S3) and (S4) is done by subroutine simplex().
4 Possible variants
4.1 Subset of R n as domain
Our experiments have shown, that the basic algorithm works even for densities
with support Since the hat h(x) has support
the rejection constant might become very big.
Subroutine 3 Construct a hat function
Input: level of triangulation for finding s, level of minimal triangulation
1: for all tuples
2: Append new cone to list of cones with en as its spanning
vectors.
3: repeat
4: for all cone C in list of cones do
5: call split() with C.
Update list of cones.
7: until level of triangulation for finding s is reached
8: for all cone C in list of cones do
9: call find() with C.
10: repeat
11: for all cone C in list of cones do
12: call split() with C.
13: Update list of cones.
14: until minimum level of triangulation is reached
15: repeat
16: for all cone C in list of cones where s unknown do =  (13) violated
17: call split() with C and list of cones.
call find() with both new cones.
19: Update list of cones.
20: until no such cone was found
21: for all cone C in list of cones do
22: Compute all parameters of C.
Total volume below hat
24: for all cone C in list of cones do
Used for O(0)-search algorithm
27: return list of cones, H tot .
Subroutine 4 cone and update list of cones
Input: cone C, list of cones
1: Find lowest indices i; j of all vectors of C.
2: Find highest index m of all vectors (of triangulation).
3:
4: Append new cone C 0 to list and copy vectors and s of C into C 0 .
5: Replace t i by t m+1 in C and replace t j by t m+1 in C 0 .
Replace det by 1
in C and C 0 .
7: return list of cones.
Subroutine 5 find a proper touching point
Input: cone C
Bracketing a minimum
1: Search for a s 1 2 DA . return failed if not successful.
2: Search for s 0 , s 2 (Use property (21)). return failed if not successful.
3: Find s using Brent's algorithm (Use (29)). return failed if not successful.
Subroutine 6 Generate a random tuple with density proportional to hat
Input: H tot , list of cones
1: Generate a uniformly [0; H tot ] distributed random variate U .
2: Find C, such that H cum
pred
C .
(C pred is the predecessor of C is the list of cones.)
3: Generate a gamma(n; fi) distributed random variate G.
uniformly distributed point in Q(G) and return tuple
4: X / call simplex() with C and G.
5: return X.
Subroutine 7 Generate a uniform distributed tuple on simplex
Input: cone C, x (location of sweeping plane)
uniformly distributed random variates U i in simplex
1: Generate iid uniform [0; 1] random variates U i ,
2: Un / 1.
3:
4: for do
5: U i / U
uniformly distributed point X in Q(x)
x
7: return X.
Pyramids. If the given domain D is a proper subset of R n (that is, we give
constraints for suppf ), the acceptance probability can be increased when we
restrict the domain of h accordingly to the domain D. (The domain is the set
of points where the density f is defined; obviously suppf ' D. Notice that we
have to provide the domain D for the algorithm but the support of f is not
known.)
Thus we replace (some) cones by pyramids. Notice that the base of such
a pyramid must be perpendicular to the direction g. Hence we first have to
choose a construction point p and then compute the height of the pyramid.
The union of these pyramids (and of the remaining cones) must cover D.
Whenever we get a random point not in the domain D we reject it. It is clear
that continued triangulation decreases the volume between D and enclosing set.
Polytopes. We only deal with the case where D is an arbitrary polytope
which are given by a set of linear inequalities.
Height of pyramid. The height is the maximum of hg; xi in C " D. Because
of our restriction to polytopes this is a linear programming problem. Using the
spanning vectors t as basis for the R n , it can be solved by means of the
simplex algorithm in at most k pivot steps (for a simple polytope), where k is
the number of constraints for D.
Marginal density and volume below hat. The marginal distribution is a
truncated gamma distribution with domain [0; u], where u is the height of the
pyramid C. Instead of (28) and (29) we find for pyramids
Z ux \Gamman \Gamma(n; fiu) (30)
and
where \Gamma(n;
R x
\Gammat dt is proportional to the incomplete gamma function
and can be computed by means of formula (3.351) in [GR65].
Computing the height u(s) is rather expensive. So it is recommended to use
instead of the exact function (31) for finding a touching point in pyramid C.
Computational experiments with the standard normal distribution have shown,
that the effect on the rejection constant is rather small (less than 5%).
4.2 Density not differentiable
For the construction of the hat function we need a tangent plane for every x 2 D.
Differentiability of the density is not really required. Thus it is sufficient to have
a subroutine that returns the normal vector of a tangent hyperplane (which is
However for densities f which are not differentiable the function in (29)
might have a nasty behavior. However notice that f must be continuous in the
interior of suppf , since log ffi f is concave.
4.3 Indicator Functions
If is the indicator function of a convex set, then we can choose an
arbitrary point in the convex set as the mode (as origin of our construction) and
set t, the center of the cone C (see (4) in x2.2). Notice that the marginal
density in (16) now reduces to h g
. None of the parameters ff and
fi depends on the choice of the touching point p. Of course we have to provide
a compact domain for the density.
Using indicator functions we can generate uniformly distributed random
variates of arbitrary convex sets.
4.4 Mode not in Origin
It is obvious that the method works, when the mode m is an arbitrary point
in R. If the mode is unknown we can use common numerical methods for
finding the maximum of f , since T (f(x)) is concave (see e.g. [Rao84]).
Notice that the exact location of the mode is not really required. The algorithm
even works if the center for the construction of the cones is not close to
the mode. Then we just get a hat with a worse rejection constant.
4.5 Add mode as construction point
Since we have only one construction point in each cone, the rejection constant
is bounded from below. Thus only a few steps to triangulate the S
sense. To get a better hat function we can use the mode m of f as an additional
construction point. The hat function is then the minimum of f(m) and the
original hat. The cone is split into two parts by a hyperplane F (b) with different
marginal densities, where b is given by T f(m). Its marginal density
is then given by
Notice that we use the same direction g for the sweep plane in both parts. We
have to compute the volume below the hat for both parts which are given by
a x
4.6 More construction points per cone
A way to improve the hat function is to use more than one (or two) construction
points. But this method has some disadvantages and it is not recommended to
use it.
overlapping region

Figure

3: Two construction points in a cone
The cones are divided in several pieces of a pyramid (see figure 3). The
lower and upper base of these pieces must be perpendicular to the corresponding
direction g. These vectors g are determined by the gradients of the transformed
density at the construction points in this pieces. Thus these g (may) differ and
hence these pieces must overlap. This increases the rejection constant. Moreover
it is not quite clear how to find such pieces. For the univariate case appropriate
methods exist (e.g. [DH94]). But in the multivariate case these are not suitable.
Also adaptive rejection sampling (introduced in [GW92]) as used in [H-or95b,
LH98] is not a really good choice. The reason is quite simple. The cones are
fixed and the construction points always are settled in the center of these cones.
Thus using adaptive rejection sampling we select the new construction points
due to a distribution which is given by the marginal density of
this marginal density is not zero at the existing construction points.
4.7 Squeezes
We can make a very simple kind of squeezes: Let x
Compute the minima of the transformed density at Q(x i ) for all i. Since ~
f is
concave these minima are at the vertices of these simplices. The squeeze s i (x)
for x
denotes the minimum of ~
f(x) in Q(x i ). The setup of these squeezes
is rather expensive and only useful, if many random points of the same distribution
must be generated.
-concave densities
A family T c of transformations that fulfil conditions (T1)-(T4) is introduced in
[H-or95a]. Let c - 0. Then we set
c (x)
It can easily be verified, that condition (T4) (i.e. volume below hat is bounded)
holds if and only if \Gamma 1
To ensure the negativity of the transformed hat we always have to choose
the mode m as construction point (see x4.5).
In [H-or95a] it was shown that if a density f is T c -concave then it is T c1 -
concave for all c 1 - c. The case log(x) is already described in
x2.6.
For the case c ! 0 the marginal density function (16) is now given by
c for x ? b
where b is given by (fi
To our knowledge no special generator
for this distribution is known. (The part for x ? b looks like a beta-prime
distribution (see [JKB95]), but ff; fi ? 0.
By assumption (fi x
\Gamman. Hence it can easily be
seen that the marginal density is T c -concave. Therefore we can use the universal
generator ([H-or95a]).
Computational Experience
5.1 A C-implementation.
A test version of the algorithm was written in C and is available via anonymous
ftp [Ley98].
It should handle the following densities
ffl f is log-concave but not constant on its support.
ffl Domain D is either R n or an arbitrary rectangle [a
ffl The mode m is arbitrary. But if D 6= suppf then m must be an interior
point of suppf not "too close" to the boundary of suppf .
We used two lists for storing the spanning vectors and the cones (with pointers
to the list of vectors). For the setup we have to store the edges (i;
computing the new vertices. This is done temporarily in a hash table, where
the first index i is used as the hash index.
The setup step is modified in the case of a rectangular domain. If the mode
is near to the boundary of D we use the nearest point on the boundary (if
possible a vertex) for the center to construct the cones. If this point is on the
boundary we easily can eliminate all those initial cones, that does not intersect
D. If this point is a vertex of the rectangle there remains only one initial cone.
For finding the mode of f we used a pattern search method by Hooke and
Jeeves [HJ61, Rao84] as implemented in [Kau63, BP66], since it could deal
with both unbounded and bounded support of f without giving explicit con-
straints. For finding the minimum of (29) we use Brent's algorithm as described
in ([FMM77]). The implementation contains some parameters to adjust these
routines to a particular density f .
For finding a cone C in subroutine sample() we used a O(0)-algorithm
with a search table. (Binary search is slower.) For generating the gamma
distributed random number G we used the algorithm in [AD82] for the case
of unbounded domain. When D is a rectangle, we used transformed density
rejection ([H-or95a]) to generate from the truncated gamma distributions. Here
it is only necessary to generate a optimal hat function for the truncated gamma
distribution with shape parameters n and 1 with domain (0; umax ), where umax
is the maximal value of height \Delta fi for all cones. The optimal touching points for
this gamma distribution are computed by means of the algorithm [DH94].
The code was written for testing various variants of the algorithm and is not
optimized for speed. Thus the data shown in the tables below give just an idea
of the performance of the algorithm.
We have tested the algorithm with various multivariate log-concave distributions
in some dimensions. All our tests have been done on a PC with a P90
processor running Linux and the GNU C compiler.
5.2 Basic version: unbounded domain, mode in origin
Random points with density proportional to hat function. The time
for the generation of random points below the hat has shown to be almost linear
in dimension n. Table 1 shows the average time for the generation of a single
point. For comparison we give the time for generation of n normal distributed
points using the Box/Muller method [BM58] (which gives a standard multi-normal
distributed point with density proportional to
For computing the hat function we only used initial cones for the standard
multinormal density.
hat function 14.6 17.1 21.3 24.9 30.2 34.6 41.5 45.7 55.6
multinormal 7.2 10.8 14.4 18.0 21.6 25.2 28.8 32.4 36.0

Table

1: average time for the generation of one random point (in -s)
Random points for the given density. The real time needed for the generation
of a random point for a given log-concave density depends on the rejection
constant and the costs for computing the density. Table 2 shows the acceptance
probabilities and the times needed for the generation of standard multinormal
distributed points. Notice that these data do not include the time for setting
up the hat function.
Setup. When find() is called after triangulation has been done, the time
needed for the computation of the hat function depends linearly on the number
number of cones
time (-s) 23.7 27.9 38.2 49.8 73.8 99.3 142 262 575
acceptance

Table

2: acceptance probability and average time for the generation of standard
multinormal distributed points
of cones. (Thus find() is the most expensive part of the setup().) Table 3
shows the situation for the multinormal distribution with density proportional
to
It demonstrates the effects of continuing
barycentric subdivision of the "oldest" edge (see x2.4) on the number of cones,
the acceptance probability, the costs for generating a random point proportional
to the hat function (i.e. without rejection) and proportional to the given density.
Furthermore it shows the total time (in ms) for the setup (i.e. for computing
the parameters of the hat function) (in ms) and the time for each cone (in -s)
(The increase for large n in the time needed for generating points below the hat
is due to effects of memory access time.)
subdivisions
cones
acceptance
hat (-s) 24.8 24.8 24.9 25.0 25.2 25.3 25.7 26.1 27.1 27.6 28.4
density (-s) 94.2 73.0 59.9 52.1 45.6 42.9 40.1 39.3 39.5 39.6 40.4
setup/cone (-s) 700 706 738 720 713 714 727 756 762 763 764

Table

3: time for computing the hat function for multinormal distribution
If we do not run find() for every cone of the triangulation but use the
method described in x3.1 we can reduce the costs for the construction of the hat
function. Table 4 gives an idea of this reduction for the multinormal distribution
with proportional to
It shows the time for
constructing the hat function subject to the number of cones for which find()
is called.
Due to

Table

4 the acceptance probability is not very bad, if we run find()
only for the initial cones. But this is not true in general. It might become extremely
poor if the level sets of the density are very "skinny". Table 5 demonstrates
the effect on the density proportional to
2.
At last table 6 demonstrates that the increase in the time for constructing
the hat function for increasing dimension n is mainly due to the increase of
find() in subdivision
cones
acceptance
setup (ms) 66.2 76.8 100.2 141.6 224.7 393.6 744
setup/cone (total) (-s)

Table

4: time for computing the hat function for multinormal distribution with
"inherited" construction points
find() in subdivision 1 2 3 4
acceptance

Table

5: acceptance probability for multinormal distribution with "inherited"
construction points
number of cones. Notice that we start with 2 n cones. Furthermore we have to
make consecutive subdivisions to shorten every edge of a simplex that
defines an initial cone. Thus the number of cones increases exponentially.
acceptance
setup/cone (-s) 540 616 714 811 927 1011 1170 1250 1421

Table

time for computing the hat function for multinormal distribution
subdivisions of the initial cones)
If the covariance matrix of the multinormal distribution is not a diagonal
matrix and the ratio of the highest and lowest eigenvalue is large, then we cannot
use initial cones only and we have to make several subdivisions of the cones.
Because of the above considerations the necessary number of cones explodes
with increasing n. Thus in this case this method cannot be used for large n.
(Suppose we have to shorten every edge of each simplex, then we have
cones if but we need 2
Tests. We ran a - 2 -test with the density proportional to exp(\Gamma
to validate the implementation. For all other densities we compared the
observed rate of acceptance to the expected acceptance probability.
Comparison with algorithm [LH98]. The code for algorithm [LH98] is
much longer (and thus contains more bugs). The setup is much slower and it
needs 11 750 -s to generate on mulitnormal distributed random point in dimension
4 (versus 38 -s in table 2 for tdrmv()).
5.3 Rectangular domain
Normal densities restricted to an arbitrary rectangle have a similar performance
as the corresponding unrestricted densities. except of the acceptance probability
which is worse since the domain of the hat h is a superset of the domain of density
f .
5.4 Quality
The quality of non-uniform random number generators using transformation
techniques is an open problem even for the univariate case (see e.g. [H-or94]
for a first approach). It depends on the underlying uniform random number
generators. The situation is more serious in the multivariate case. Notice that
this new algorithm requires more than n+2 uniform random numbers for every
random point. We cannot give an answer to this problem here, but it should be
clear that e.g. RANDU (formerly part of IBM's Scientific Subroutine Package,
and now famous for its devasting defect in three dimensions: its consecuting
points just fifteen parallel planes; see e.g. [LW97]) may
result in a generator of poor quality.
5.5 Some Examples
We have tested our algorithm in dimensions
proportional to
where a i ? 0. The domain was R n and some rectangles. We also used densities
proportional to f i (U x U is an orthonormal transformation and b a
vector, to test distributions with non-diagonal correlation matrix and arbitrary
mode.
The algorithm works well for densities f 3 , f 4 and f 5 both with
and D being a rectangle enclosing the support of f i . Although some of these
densities are not C 1 , the find() routine works. Problems arise if the level sets
of the density have "corners", i.e. the g is unstable when we vary the touching
a little bit. Then there are somes (that contains these "corners") with huge
volumeHC and further triangulation is necessary. If dimension is high (n & 5)
too many cones might be necessary. The optimization algorithm for finding the
mode fails if we use a starting point outside the support of f 5 .
The code has some parameters for adjusting the algorithm to the given
density. For example, it requires some testing to get the optimal number of
cones and the optimal level of subdivisions for calling find().
5.6 R'esum'e
The presented algorithm is a suitable method for sampling from log-concave
(and T -concave) distributions. The algorithm works well for all tested log-concave
densities if dimension is low (n . 5) or if correlation is not too high.
Restrictions of these densities to compact polytopes are possible. The setup
time is small for small dimension but increases exponentially in n. The speed
for generating random points is quite fast even for n - 6. Due to the large
amount of cones for high dimension the program requires a lot of computer
memory (typically 2-10 MB).
Although the developed algorithm is not a real black box method it is easily
adjustable for a large class of log-concave densities. Examples for which
the algorithm works are the multivariate normal distribution and the multi-variate
student distribution (with transformation T arbitrary
mean vector and variance matrix conditioned to an arbitrary compact polytope.
However for higher dimensions the ratio of highest and lowest eigenvalue of the
covariance matrix should not be "too big".

Acknowledgments

The author wishes to note his appreciation for help rendered by J-org Lenneis.
He has given lots of hints for the implementation of the algorithm. The author
also thanks Gerhard Derflinger and Wolfgang H-ormann for helpful conversations
and their interest in his work.



--R

Generating gamma variates by a modified rejection technique.

Box and M.
Remark on algorithm 178.
Principles of Random Variate Generation.

Random variate generation for multivariate densities.
The optimal selection of hat functions for rejection algorithms.
Random variable generation using concavity properties of transformed densities.
Computer methods for mathematical computations.
Table of Integrals.
Convex Poytopes.
Adaptive rejection sampling for gibbs sam- pling
Universal generators for correlation induction.
"Direct search"
The quality of non-uniform random numbers
A rejection technique for sampling from T
A universal generator for bivariate log-concave distri- butions
Continuous Univariate Distributions
Multivariate Statistical Simulation.

Polytope Volume

A sweep plane algorithm for generating random tuples.
Inversive and linear congruential pseudorandom number generators in empirical tests.

Theory and Applications.
On computer generation of random vectors by transformations of uniformly distributed vectors.
The Computation of Fixed Points and Applications
Improving the convergence of fixed-point algorithms
Efficient generation of random variates via the ratio-of-uniforms method
Lectures on Polytopes
--TR
Multivariate statistical simulation
On computer generation of random vectors by transformation of uniformly distributed vectors
A rejection technique for sampling from <italic>T</italic>-concave distributions
Inversive and linear congruential pseudorandom number generators in empirical tests
Random variate generation for multivariate unimodal densities
A sweep-plane algorithm for generating random tuples in simple polytopes
`` Direct Search'''' Solution of Numerical and Statistical Problems
Generating gamma variates by a modified rejection technique
Remark on algorithm 178 [E4] direct search
Algorithm 178: direct search
Computer Methods for Mathematical Computations

--CTR
G. Leobacher , F. Pillichshammer, A method for approximate inversion of the hyperbolic CDF, Computing, v.69 n.4, p.291-303, December 2002
Wolfgang Hrmann, Algorithm 802: an automatic generator for bivariate log-concave distributions, ACM Transactions on Mathematical Software (TOMS), v.26 n.1, p.201-219, March 2000
W. Hrmann , J. Leydold, Random-number and random-variate generation: automatic random variate generation for simulation input, Proceedings of the 32nd conference on Winter simulation, December 10-13, 2000, Orlando, Florida
Seyed Taghi Akhavan Niaki , Babak Abbasi, Norta and neural networks based method to generate random vectors with arbitrary marginal distributions and correlation matrix, Proceedings of the 17th IASTED international conference on Modelling and simulation, p.234-239, May 24-26, 2006, Montreal, Canada
sampling with the ratio-of-uniforms method, ACM Transactions on Mathematical Software (TOMS), v.26 n.1, p.78-98, March 2000
rejection method;multivariate log-concave distributions