- Email: [email protected]
Generating generalized inverse Gaussian random variates by fast inversion
Computational Statistics and Data Analysis 55 (2011) 213–217
Contents lists available at ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Generating generalized inverse Gaussian random variates by fast inversionI Josef Leydold a,∗ , Wolfgang Hörmann b a
Department of Statistics and Mathematics, WU (Vienna University of Economics and Business), Augasse 2–6, A-1090 Wien, Austria
b
Department of Industrial Engineering, Boğaziçi University, 34342 Bebek-İstanbul, Turkey
article
info
Article history: Received 18 November 2009 Received in revised form 12 July 2010 Accepted 13 July 2010 Available online 18 July 2010 Keywords: Generalized inverse Gaussian distribution Random variate generation Numerical inversion
abstract The inversion method for generating non-uniformly distributed random variates is a crucial part in many applications of Monte Carlo techniques, e.g., when low discrepancy sequences or copula based models are used. Unfortunately, closed form expressions of quantile functions of important distributions are often not available. The (generalized) inverse Gaussian distribution is a prominent example. It is shown that algorithms that are based on polynomial approximation are well suited for this distribution. Their precision is close to machine precision and they are much faster than root finding methods like the bisection method that has been recently proposed. © 2010 Elsevier B.V. All rights reserved.
1. Introduction There are two recent trends in the framework of financial engineering. On the one hand, the development of more realistic models leads to an increasing number of applications of less frequently used distributions. On the other hand, the complexity of these new models requires more sophisticated numerical algorithms beyond naïve Monte Carlo methods. Among these copula based methods, the usage of low discrepancy sequences (the so-called quasi-Monte Carlo methods) is of great importance. However, these methods require the inversion method, i.e., the monotone one-to-one transformation of uniform (0, 1) pseudo-random numbers into random variates from the requested distribution. Unfortunately, the distributions in these new models are often nasty to handle. In particular, computing the inverse of the cumulative distribution function (CDF) is extremely expensive. This is the case for the generalized inverse Gaussian (GIG) distribution. Its density is given by
1 χ (ψ/χ )θ/2 θ−1 x exp − + ψx , fgig (x; θ , ψ, χ ) = 2Kθ (ψχ ) 2 x 0,
x > 0,
(1)
x ≤ 0,
where Kθ (·) denotes the modified Bessel function of the third kind with index θ (Jørgensen, 1982; Abramowitz, 1972). The domain of variation for the parameters is given by
θ ∈ R,
( {(ψ, χ ): ψ > 0, χ ≥ 0} (ψ, χ ) ∈ {(ψ, χ ): ψ > 0, χ ≥ 0} {(ψ, χ ): ψ ≥ 0, χ ≥ 0}
if θ > 0, if θ = 0, if θ < 0.
I Examples and programs appear as annex in the electronic version of this paper.
∗
Corresponding author. Tel.: +43 1 31336 4695; fax: +43 1 31336 774. E-mail addresses: [email protected] (J. Leydold), [email protected] (W. Hörmann).
0167-9473/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2010.07.011
214
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
The inverse Gaussian distribution (IG, also called the Wald distribution) is a special case of (1) with θ = − 21 . Its density is often written as
r fig (x; µ, λ) =
λ λ(x − µ)2 exp − , 2π x3 2µ2 x
x; µ, λ > 0.
(2)
Its CDF can be written as
r r ! ! λ x 2λ λ x −1 + exp Φ − +1 , Fig (x; µ, λ) = Φ x µ µ x µ
(3)
where Φ (·) denotes the standard Gaussian distribution function. For such distributions the inverse CDF is not known explicitly. Lai (2009) thus proposes to solve the nonlinear equation U = F (X ) by means of the bisection method (he specifically uses representation to implement the bisection method). He reports that the Newton–Raphson method does not always converge in his tests (with different parameter settings). The bisection method converges slowly, however, and is thus usually considered as a last resort. For the GIG distribution, the situation is even worse since a numeric quadrature routine is needed. For the case of fixed parameters, we propose two much faster methods. These methods were presented in two papers (Hörmann and Leydold, 2003; Derflinger et al., 2010), and are based on polynomial approximations to the inverse CDF of the (generalized) inverse Gaussian distribution. In the remaining part of this note, we shortly describe these algorithms and report the computational results. 2. Inversion by polynomial interpolation The basic idea of this approach is quite simple: Evaluate F , the CDF, at a couple of points xi and interpolate the nodes
(ui = F (xi ), xi ) by means of polynomials. This is done in the setup part of the algorithm where all required coefficients are computed and stored in a table. Notice that there is no need to compute F −1 (ui ) during the setup. This table is then used in the sampling part to compute the approximate inverse CDF for a given U. There exist several variants for this task. In Hörmann and Leydold (2003) we argued in favor of (cubic) Hermite interpolation for the following reasons:
• The approximation error can be bounded for smooth CDFs. • It allows one to estimate the approximation error ‘‘on the fly’’ during the setup. The maximal tolerated error can thus be controlled by the user.
• It is possible to adjust the interpolant stepwise by adding additional points without recomputing all polynomials. We thus start with a rough approximation with only a few nodes and then recursively add more nodes wherever improvements are needed. • The coefficients of the polynomials need to satisfy very simple conditions to ensure monotonicity of the interpolant. Imai and Tan (2009) report that this method performs very well for the generalized hyperbolic distribution and the normal–inverse Gaussian distribution. They use numerical integration to compute the CDF, however, and report a slow setup time. In Derflinger et al. (2010) we suggest a new algorithm that combines Newton’s interpolation formula with adaptive Gauss–Lobatto integration to evaluate the differences of the CDF. The probability density function (PDF) of the target distribution is thus the only input required. The resulting setup is much faster than using the algorithm of Hörmann and Leydold (2003) together with a packed quadrature routine like QUADPACK (Piessens et al., 1983) to compute the CDF. Moreover, we observed that it is numerically more robust when the maximal tolerated approximation error is close to machine precision. A main concern of any numerical inversion algorithm must be the control of the approximation error, i.e., the deviation of the approximate inverse CDF, Fa−1 , from the exact function F −1 . We are convinced that the u-error defined by
εu (u) = |u − F (Fa−1 (u))|
(4)
is well suited for this task. In particular, it can be computed during the setup and it can be interpreted with respect to the resolution of the underlying uniform pseudo-random number generator or low discrepancy set (see Derflinger et al., 2010, for details). In fact, goodness-of-fit tests like the Kolmogorov–Smirnov test or the χ 2 test look at that deviation. In what follows, we call the maximal tolerated u-error the u-resolution of the algorithm. The number of required nodes (and thus table size and setup time) depends on the requested u-resolution and on the target distribution. The marginal generation time is very fast and almost independent from the target distribution (the runtime is only influenced by the table size due to cache effects).
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
215
3. Computational experiments We have worked out all required details and implemented these two algorithms in our C library UNU.RAN (Leydold and Hörmann, 2010a). We used the names HINV for the Hermite interpolation algorithm of Hörmann and Leydold (2003) and PINV for the Newton interpolation with Gauss–Lobatto integration algorithm of Derflinger et al. (2010). The UNU.RAN library is also accessible from the R statistical programming language (R Development Core Team, 2009) by using our package Runuran (Leydold and Hörmann, 2010b). The following example shows how method PINV can be used in R to draw a sample from an IG distribution. Notice that the setup and sampling part are separated. The tables for the interpolant are stored in object gen. Thus there is no necessity to rerun the setup every time.
> > > > > > >
## load ‘Runuran’ library library(Runuran) ## define PDF for inverse Gaussian ig.pdf <- function(x,mu,lambda) { sqrt(1/x^{3}) * exp(-(lambda*(x-mu)^{2})/(2*mu^{2}*x)) } ## run setup (create generator object for mu=3 and lambda=2) gen <- pinv.new(pdf=ig.pdf, lb=0, ub=Inf, center=1, uresolution=1.e-10, mu=3, lambda=2) > ## draw a sample of size 5 > ur(gen,5) [1] 0.6151017 3.5961721 0.7396323 0.9170969 1.7725783 > ## compute (approximate) inverse CDF for u=0.9 > uq(gen,0.9) [1] 6.84101 The argument uresolution is optional; if omitted the u-resolution is initialized to 10−10 . The argument center is a point in the central part of the distribution and is merely a hint to the algorithm where the most relevant part of the domain is located. (The center is denoted by xc in Algorithm NINIGL in Derflinger et al. (2010).) The mode or a point near the mode is a good choice. More details about the above code can be found in the Runuran user manual. Remark 1. Our algorithms are designed as black-box algorithms, i.e., the user can provide a function that evaluates the PDF and a maximal tolerated approximation error. Our example could thus be modified to suit any distribution. (Of course it does not work for all distributions, but we observed no problems for densities that are bounded, smooth, and have moderate to low tails.) Coding the PDF can be avoided as the package already has built-in distributions. Thus udig(mu=3, an object that contains all required information (including the PDF) about the IG distribution.
lambda=2) creates
> ## run setup (create generator object for IG with mu=3 and lambda=2) > gen <- pinvd.new( udig(mu=3, lambda=2), uresolution=1.e-12 ) > ## draw a sample of size 5 > ur(gen,5) [1] 8.2701337 0.9978058 1.1816979 0.4788929 2.0146778 It is not astonishing that this method is also suitable for sampling from the generalized inverse Gaussian distribution with the density given in (1).
> gen <- pinvd.new( udgig(theta=-1, psi=3, chi=2) ) > ur(gen,5) [1] 0.1987231 0.8583419 1.4482379 0.2322391 1.6078295 Our extensive tests with the IG distribution and CDF (3) (we considered many different parameters) show that the observed u-error never exceeds a requested u-resolution of 10−10 or greater. For a u-resolution of 10−12 , we came across a few parameter settings where the maximal observed u-error was slightly too large, but never larger than 1.16 · 10−12 . We also measured both the setup time and the marginal generation time for our methods PINV and HINV as well as for the bisection method BIS. In addition we added a modified regula falsi algorithm RF (Neumaier, 2001) that uses bisection when convergence is too slow. For both RF and BIS we used the 10th and 90th percentiles for the starting interval as well as a table of size 100 (a table of size 1000 did not show much improvement). All methods are available in our UNU.RAN library. Table 1 summarizes our timing results. The reported timings are relative to sampling one exponential distributed random variate using inversion (i.e., -log(1-U)). As can easily be seen, the methods based on polynomial interpolation have very fast marginal generation times (faster than generating an exponential random variate) but require some setup. This is worthwhile for moderate (or larger) sample sizes; break-even points between PINV and the two root finding algorithms are listed in Table 2. Regula falsi is faster only when samples smaller than that break-even points are required. We could not see
216
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
Table 1 Marginal generation time and setup time (in parenthesis) for some IG(µ, λ) and GIG(θ, ψ, χ ) distributions and u-resolutions εu . The timings are relative to the sampling time for one exponential distributed random variate using inversion (i.e., -log(1-U)). PINV-5: Newton interpolation of order 5. HINV-3: cubic Hermite interpolation. RF-2 & RF-100: modified regula falsi with a table of starting points of sizes 2 and 100. BIS-2 & BIS-100 bisection method with a table of starting points of sizes 2 and 100. The CDF for GIG is computed by numerical integration (the runtimes are indicated in italics). Time unit is 0.105 µs., sample size 105 ; Intel Core Duo 2.0 GHz, Linux 2.6.26, GCC 4.3.4. Method
εu
PINV-5
10−8 10−10 10−12
0.53 (12 651) 0.54 (20 002) 0.53 (33 859)
0.54 (8 354) 0.53 (15 677) 0.53 (29 311)
0.53 (15 935) 0.53 (23 118) 0.53 (38 110)
0.53 0.53 0.53
HINV-3
10−8 10−10 10−12
0.65 (4 539) 0.66 (12 893) 0.74 (41 409)
0.66 (4 195) 0.66 (13 012) 0.74 (41 770)
0.66 (4 443) 0.68 (12 886) 0.75 (41 282)
0.56 (332 021) 0.56 (717 424) 0.57 (1856 165)
RF-2
10−8 10−10 10−12
25.53 27.63 29.05
RF-100
10−8 10−10 10−12
16.96 (2 182) 18.45 (2 339) 20.16 (2 557)
BIS-2
10−8 10−10 10−12
87.62 111.96 136.39
BIS-100
10−8 10−10 10−12
69.64 (2 187) 94.05 (2 325) 118.45 (2 548)
IG (0.45, 0.3)
(131) (130) (137)
(117) (106) (115)
IG (0.5, 18)
20.92 22.31 23.61
(124) (129) (138)
13.93 (1 824) 15.76 (1 963) 16.97 (2 136) 75.38 96.53 117.57
(126) (139) (131)
59.60 (1 825) 80.71 (1 985) 101.91 (2 150)
IG (4.3, 0.6)
23.25 24.96 26.16
GIG (−5, 2, 3) (11 374) (19 366) (37 043)
(95) (100) (85)
923 1008 1057
(8 854) (8 585) (8 225)
14.83 (1 810) 15.91 (2 041) 17.80 (2 119)
677 750 800
(84 552) (93 848) (97 676)
(88) (74) (130)
3570 4595 5597
(9 733) (8 440) (10 081)
55.72 (1 817) 75.00 (2 027) 94.09 (2 120)
2894 3914 4922
(84 052) (90 329) (96 453)
69.52 88.65 107.84
Table 2 Break-even points between PINV and the methods RF and BIS, respectively. Method
εu
IG (0.45, 0.3)
IG (0.5, 18)
IG (4.3, 0.6)
GIG (−5, 2, 3)
RF-2
10−8 10−10 10−12
520 725 1218
394 705 1246
687 928 1462
3 12 27
BIS-2
10−8 10−10 10−12
150 178 258
108 159 244
228 255 352
3 4 7
any benefit from using the bisection method as it was always slower than regula falsi. We also tried numerical integration (using the GNU Scientific library Galassi et al., 2009) to obtain the CDF for the GIG distribution. However, the runtimes are very slow (see Table 1 for details). We also tested Algorithm PINV with the inverse Gaussian distribution with small values for mean µ or shape parameter λ. In particular, we ran experiments for the IG distribution for µ = 1 and λ = 10−1 , 10−2 , . . . , 10−7 as well as with parameters λ = 1 and µ = 10−1 , . . . , 10−7 . Down to a desired u-resolution of 10−12 the algorithm did not show any problems. In one case we had to adapt the domain of the distribution. Remark 2. Our computational experiments (as well as those of Imai and Tan, 2009) make it clear that the problems of HINV (‘‘HL method’’) reported by Lai (2009) in his last two paragraphs can only be attributed to a problem in his implementation. Indeed parameter δ of the ‘‘HL method’’ is used for the initial rough approximation and has hardly any influence on the accuracy of the resulting interpolant (once the setup is terminated). Also compare our short, simple, and slow R implementation of HINV contained as file HINV_Rcode.R in the electronic supplement of this paper which can be found online at doi:10.1016/j.csda.2010.07.011. That implementation exactly follows the algorithm description in Hörmann and Leydold (2003). The examples within the same file demonstrate that for the normal, exponential, gamma, and inverse Gaussian distributions, HINV reaches the desired u-resolution without problems. 4. Conclusion Our algorithms are well suited for generating random variates from the (generalized) inverse Gaussian distribution by numerical inversion. In all our experiments the maximal observed u-error was always within acceptable deviation from the requested u-resolution. After a moderate setup, our algorithms are faster than the standard method for generating exponential variates and much faster than numerical root finding algorithms such as the bisection method that has been suggested by Lai (2009). The speed-up factors we observed are between 50 and 220 for sample sizes 105 or larger. Our algorithms are available as C library UNU.RAN and as R package Runuran which can both be downloaded from our web site.
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
217
Acknowledgement The authors are grateful to Carsten Botts for his useful and valuable comments. Appendix. Supplementary data Supplementary material related to this article can be found online at doi:10.1016/j.csda.2010.07.011. References Abramowitz, M., Stegun, I.A. (Eds.), 1972. Handbook of Mathematical Functions, 9th ed.. Dover, New York, ISBN: 0486-612-724. Derflinger, G., Hörmann, W., Leydold, J., 2010. Random variate generation by numerical inversion when only the density is known. ACM Trans. Model. Comput. Simul 20 (4), #18. Galassi, M., et al. GNU Scientific Library. Version 1.13. 2009. http://www.gnu.org/software/gsl/. Hörmann, W., Leydold, J., 2003. Continuous random variate generation by fast numerical inversion. ACM Trans. Model. Comput. Simul. 13 (4), 347–362. Imai, J., Tan, K.S., 2009. an accelerating quasi-Monte Carlo method for option pricing under the generalized hyperbolic Lévy process. SIAM J. Sci. Comput. 31 (3), 2282–2302. doi:10.1137/080727713. Jørgensen, B., 1982. Statistical Properties of the Generalized Inverse Gaussian Distribution. In: Lecture Notes in Statistics, vol. 9. Springer-Verlag, New York, Berlin. Lai, Y., 2009. Generating inverse Gaussian random variates by approximation. Comput. Statist. Data Anal. 53 (10), 3553–3559. doi:10.1016/j.csda.2009.03.007. Leydold, J., Hörmann, W., 2010a. UNU.RAN—a library for non-uniform universal random variate generation. Version 1.6. Department of Statistics and Mathematics, WU Wien, A-1090 Wien, Austria. http://statmath.wu-wien.ac.at/unuran/. Leydold, J., Hörmann, W., 2010b. Runuran—R interface to the UNU.RAN random variate generators. Version 0.13. Department of Statistics and Mathematics, WU Wien, A-1090 Wien, Austria. http://cran.r-project.org/. Neumaier, A., 2001. Introduction to Numerical Analysis. Cambridge University Press, Cambridge. Piessens, R., de Doncker-Kapenga, E., Überhuber, C.W., Kahaner, D.K., 1983. QUADPACK. A Subroutine Package for Automatic Integration. In: Springer Series in Computational Mathematics, vol. 1. Springer-Verlag, Berlin. R Development Core Team. 2009. R: a language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. http://www.R-project.org. ISBN: 3-900051-07-0.
Contents lists available at ScienceDirect
Computational Statistics and Data Analysis journal homepage: www.elsevier.com/locate/csda
Generating generalized inverse Gaussian random variates by fast inversionI Josef Leydold a,∗ , Wolfgang Hörmann b a
Department of Statistics and Mathematics, WU (Vienna University of Economics and Business), Augasse 2–6, A-1090 Wien, Austria
b
Department of Industrial Engineering, Boğaziçi University, 34342 Bebek-İstanbul, Turkey
article
info
Article history: Received 18 November 2009 Received in revised form 12 July 2010 Accepted 13 July 2010 Available online 18 July 2010 Keywords: Generalized inverse Gaussian distribution Random variate generation Numerical inversion
abstract The inversion method for generating non-uniformly distributed random variates is a crucial part in many applications of Monte Carlo techniques, e.g., when low discrepancy sequences or copula based models are used. Unfortunately, closed form expressions of quantile functions of important distributions are often not available. The (generalized) inverse Gaussian distribution is a prominent example. It is shown that algorithms that are based on polynomial approximation are well suited for this distribution. Their precision is close to machine precision and they are much faster than root finding methods like the bisection method that has been recently proposed. © 2010 Elsevier B.V. All rights reserved.
1. Introduction There are two recent trends in the framework of financial engineering. On the one hand, the development of more realistic models leads to an increasing number of applications of less frequently used distributions. On the other hand, the complexity of these new models requires more sophisticated numerical algorithms beyond naïve Monte Carlo methods. Among these copula based methods, the usage of low discrepancy sequences (the so-called quasi-Monte Carlo methods) is of great importance. However, these methods require the inversion method, i.e., the monotone one-to-one transformation of uniform (0, 1) pseudo-random numbers into random variates from the requested distribution. Unfortunately, the distributions in these new models are often nasty to handle. In particular, computing the inverse of the cumulative distribution function (CDF) is extremely expensive. This is the case for the generalized inverse Gaussian (GIG) distribution. Its density is given by
1 χ (ψ/χ )θ/2 θ−1 x exp − + ψx , fgig (x; θ , ψ, χ ) = 2Kθ (ψχ ) 2 x 0,
x > 0,
(1)
x ≤ 0,
where Kθ (·) denotes the modified Bessel function of the third kind with index θ (Jørgensen, 1982; Abramowitz, 1972). The domain of variation for the parameters is given by
θ ∈ R,
( {(ψ, χ ): ψ > 0, χ ≥ 0} (ψ, χ ) ∈ {(ψ, χ ): ψ > 0, χ ≥ 0} {(ψ, χ ): ψ ≥ 0, χ ≥ 0}
if θ > 0, if θ = 0, if θ < 0.
I Examples and programs appear as annex in the electronic version of this paper.
∗
Corresponding author. Tel.: +43 1 31336 4695; fax: +43 1 31336 774. E-mail addresses: [email protected] (J. Leydold), [email protected] (W. Hörmann).
0167-9473/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2010.07.011
214
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
The inverse Gaussian distribution (IG, also called the Wald distribution) is a special case of (1) with θ = − 21 . Its density is often written as
r fig (x; µ, λ) =
λ λ(x − µ)2 exp − , 2π x3 2µ2 x
x; µ, λ > 0.
(2)
Its CDF can be written as
r r ! ! λ x 2λ λ x −1 + exp Φ − +1 , Fig (x; µ, λ) = Φ x µ µ x µ
(3)
where Φ (·) denotes the standard Gaussian distribution function. For such distributions the inverse CDF is not known explicitly. Lai (2009) thus proposes to solve the nonlinear equation U = F (X ) by means of the bisection method (he specifically uses representation to implement the bisection method). He reports that the Newton–Raphson method does not always converge in his tests (with different parameter settings). The bisection method converges slowly, however, and is thus usually considered as a last resort. For the GIG distribution, the situation is even worse since a numeric quadrature routine is needed. For the case of fixed parameters, we propose two much faster methods. These methods were presented in two papers (Hörmann and Leydold, 2003; Derflinger et al., 2010), and are based on polynomial approximations to the inverse CDF of the (generalized) inverse Gaussian distribution. In the remaining part of this note, we shortly describe these algorithms and report the computational results. 2. Inversion by polynomial interpolation The basic idea of this approach is quite simple: Evaluate F , the CDF, at a couple of points xi and interpolate the nodes
(ui = F (xi ), xi ) by means of polynomials. This is done in the setup part of the algorithm where all required coefficients are computed and stored in a table. Notice that there is no need to compute F −1 (ui ) during the setup. This table is then used in the sampling part to compute the approximate inverse CDF for a given U. There exist several variants for this task. In Hörmann and Leydold (2003) we argued in favor of (cubic) Hermite interpolation for the following reasons:
• The approximation error can be bounded for smooth CDFs. • It allows one to estimate the approximation error ‘‘on the fly’’ during the setup. The maximal tolerated error can thus be controlled by the user.
• It is possible to adjust the interpolant stepwise by adding additional points without recomputing all polynomials. We thus start with a rough approximation with only a few nodes and then recursively add more nodes wherever improvements are needed. • The coefficients of the polynomials need to satisfy very simple conditions to ensure monotonicity of the interpolant. Imai and Tan (2009) report that this method performs very well for the generalized hyperbolic distribution and the normal–inverse Gaussian distribution. They use numerical integration to compute the CDF, however, and report a slow setup time. In Derflinger et al. (2010) we suggest a new algorithm that combines Newton’s interpolation formula with adaptive Gauss–Lobatto integration to evaluate the differences of the CDF. The probability density function (PDF) of the target distribution is thus the only input required. The resulting setup is much faster than using the algorithm of Hörmann and Leydold (2003) together with a packed quadrature routine like QUADPACK (Piessens et al., 1983) to compute the CDF. Moreover, we observed that it is numerically more robust when the maximal tolerated approximation error is close to machine precision. A main concern of any numerical inversion algorithm must be the control of the approximation error, i.e., the deviation of the approximate inverse CDF, Fa−1 , from the exact function F −1 . We are convinced that the u-error defined by
εu (u) = |u − F (Fa−1 (u))|
(4)
is well suited for this task. In particular, it can be computed during the setup and it can be interpreted with respect to the resolution of the underlying uniform pseudo-random number generator or low discrepancy set (see Derflinger et al., 2010, for details). In fact, goodness-of-fit tests like the Kolmogorov–Smirnov test or the χ 2 test look at that deviation. In what follows, we call the maximal tolerated u-error the u-resolution of the algorithm. The number of required nodes (and thus table size and setup time) depends on the requested u-resolution and on the target distribution. The marginal generation time is very fast and almost independent from the target distribution (the runtime is only influenced by the table size due to cache effects).
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
215
3. Computational experiments We have worked out all required details and implemented these two algorithms in our C library UNU.RAN (Leydold and Hörmann, 2010a). We used the names HINV for the Hermite interpolation algorithm of Hörmann and Leydold (2003) and PINV for the Newton interpolation with Gauss–Lobatto integration algorithm of Derflinger et al. (2010). The UNU.RAN library is also accessible from the R statistical programming language (R Development Core Team, 2009) by using our package Runuran (Leydold and Hörmann, 2010b). The following example shows how method PINV can be used in R to draw a sample from an IG distribution. Notice that the setup and sampling part are separated. The tables for the interpolant are stored in object gen. Thus there is no necessity to rerun the setup every time.
> > > > > > >
## load ‘Runuran’ library library(Runuran) ## define PDF for inverse Gaussian ig.pdf <- function(x,mu,lambda) { sqrt(1/x^{3}) * exp(-(lambda*(x-mu)^{2})/(2*mu^{2}*x)) } ## run setup (create generator object for mu=3 and lambda=2) gen <- pinv.new(pdf=ig.pdf, lb=0, ub=Inf, center=1, uresolution=1.e-10, mu=3, lambda=2) > ## draw a sample of size 5 > ur(gen,5) [1] 0.6151017 3.5961721 0.7396323 0.9170969 1.7725783 > ## compute (approximate) inverse CDF for u=0.9 > uq(gen,0.9) [1] 6.84101 The argument uresolution is optional; if omitted the u-resolution is initialized to 10−10 . The argument center is a point in the central part of the distribution and is merely a hint to the algorithm where the most relevant part of the domain is located. (The center is denoted by xc in Algorithm NINIGL in Derflinger et al. (2010).) The mode or a point near the mode is a good choice. More details about the above code can be found in the Runuran user manual. Remark 1. Our algorithms are designed as black-box algorithms, i.e., the user can provide a function that evaluates the PDF and a maximal tolerated approximation error. Our example could thus be modified to suit any distribution. (Of course it does not work for all distributions, but we observed no problems for densities that are bounded, smooth, and have moderate to low tails.) Coding the PDF can be avoided as the package already has built-in distributions. Thus udig(mu=3, an object that contains all required information (including the PDF) about the IG distribution.
lambda=2) creates
> ## run setup (create generator object for IG with mu=3 and lambda=2) > gen <- pinvd.new( udig(mu=3, lambda=2), uresolution=1.e-12 ) > ## draw a sample of size 5 > ur(gen,5) [1] 8.2701337 0.9978058 1.1816979 0.4788929 2.0146778 It is not astonishing that this method is also suitable for sampling from the generalized inverse Gaussian distribution with the density given in (1).
> gen <- pinvd.new( udgig(theta=-1, psi=3, chi=2) ) > ur(gen,5) [1] 0.1987231 0.8583419 1.4482379 0.2322391 1.6078295 Our extensive tests with the IG distribution and CDF (3) (we considered many different parameters) show that the observed u-error never exceeds a requested u-resolution of 10−10 or greater. For a u-resolution of 10−12 , we came across a few parameter settings where the maximal observed u-error was slightly too large, but never larger than 1.16 · 10−12 . We also measured both the setup time and the marginal generation time for our methods PINV and HINV as well as for the bisection method BIS. In addition we added a modified regula falsi algorithm RF (Neumaier, 2001) that uses bisection when convergence is too slow. For both RF and BIS we used the 10th and 90th percentiles for the starting interval as well as a table of size 100 (a table of size 1000 did not show much improvement). All methods are available in our UNU.RAN library. Table 1 summarizes our timing results. The reported timings are relative to sampling one exponential distributed random variate using inversion (i.e., -log(1-U)). As can easily be seen, the methods based on polynomial interpolation have very fast marginal generation times (faster than generating an exponential random variate) but require some setup. This is worthwhile for moderate (or larger) sample sizes; break-even points between PINV and the two root finding algorithms are listed in Table 2. Regula falsi is faster only when samples smaller than that break-even points are required. We could not see
216
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
Table 1 Marginal generation time and setup time (in parenthesis) for some IG(µ, λ) and GIG(θ, ψ, χ ) distributions and u-resolutions εu . The timings are relative to the sampling time for one exponential distributed random variate using inversion (i.e., -log(1-U)). PINV-5: Newton interpolation of order 5. HINV-3: cubic Hermite interpolation. RF-2 & RF-100: modified regula falsi with a table of starting points of sizes 2 and 100. BIS-2 & BIS-100 bisection method with a table of starting points of sizes 2 and 100. The CDF for GIG is computed by numerical integration (the runtimes are indicated in italics). Time unit is 0.105 µs., sample size 105 ; Intel Core Duo 2.0 GHz, Linux 2.6.26, GCC 4.3.4. Method
εu
PINV-5
10−8 10−10 10−12
0.53 (12 651) 0.54 (20 002) 0.53 (33 859)
0.54 (8 354) 0.53 (15 677) 0.53 (29 311)
0.53 (15 935) 0.53 (23 118) 0.53 (38 110)
0.53 0.53 0.53
HINV-3
10−8 10−10 10−12
0.65 (4 539) 0.66 (12 893) 0.74 (41 409)
0.66 (4 195) 0.66 (13 012) 0.74 (41 770)
0.66 (4 443) 0.68 (12 886) 0.75 (41 282)
0.56 (332 021) 0.56 (717 424) 0.57 (1856 165)
RF-2
10−8 10−10 10−12
25.53 27.63 29.05
RF-100
10−8 10−10 10−12
16.96 (2 182) 18.45 (2 339) 20.16 (2 557)
BIS-2
10−8 10−10 10−12
87.62 111.96 136.39
BIS-100
10−8 10−10 10−12
69.64 (2 187) 94.05 (2 325) 118.45 (2 548)
IG (0.45, 0.3)
(131) (130) (137)
(117) (106) (115)
IG (0.5, 18)
20.92 22.31 23.61
(124) (129) (138)
13.93 (1 824) 15.76 (1 963) 16.97 (2 136) 75.38 96.53 117.57
(126) (139) (131)
59.60 (1 825) 80.71 (1 985) 101.91 (2 150)
IG (4.3, 0.6)
23.25 24.96 26.16
GIG (−5, 2, 3) (11 374) (19 366) (37 043)
(95) (100) (85)
923 1008 1057
(8 854) (8 585) (8 225)
14.83 (1 810) 15.91 (2 041) 17.80 (2 119)
677 750 800
(84 552) (93 848) (97 676)
(88) (74) (130)
3570 4595 5597
(9 733) (8 440) (10 081)
55.72 (1 817) 75.00 (2 027) 94.09 (2 120)
2894 3914 4922
(84 052) (90 329) (96 453)
69.52 88.65 107.84
Table 2 Break-even points between PINV and the methods RF and BIS, respectively. Method
εu
IG (0.45, 0.3)
IG (0.5, 18)
IG (4.3, 0.6)
GIG (−5, 2, 3)
RF-2
10−8 10−10 10−12
520 725 1218
394 705 1246
687 928 1462
3 12 27
BIS-2
10−8 10−10 10−12
150 178 258
108 159 244
228 255 352
3 4 7
any benefit from using the bisection method as it was always slower than regula falsi. We also tried numerical integration (using the GNU Scientific library Galassi et al., 2009) to obtain the CDF for the GIG distribution. However, the runtimes are very slow (see Table 1 for details). We also tested Algorithm PINV with the inverse Gaussian distribution with small values for mean µ or shape parameter λ. In particular, we ran experiments for the IG distribution for µ = 1 and λ = 10−1 , 10−2 , . . . , 10−7 as well as with parameters λ = 1 and µ = 10−1 , . . . , 10−7 . Down to a desired u-resolution of 10−12 the algorithm did not show any problems. In one case we had to adapt the domain of the distribution. Remark 2. Our computational experiments (as well as those of Imai and Tan, 2009) make it clear that the problems of HINV (‘‘HL method’’) reported by Lai (2009) in his last two paragraphs can only be attributed to a problem in his implementation. Indeed parameter δ of the ‘‘HL method’’ is used for the initial rough approximation and has hardly any influence on the accuracy of the resulting interpolant (once the setup is terminated). Also compare our short, simple, and slow R implementation of HINV contained as file HINV_Rcode.R in the electronic supplement of this paper which can be found online at doi:10.1016/j.csda.2010.07.011. That implementation exactly follows the algorithm description in Hörmann and Leydold (2003). The examples within the same file demonstrate that for the normal, exponential, gamma, and inverse Gaussian distributions, HINV reaches the desired u-resolution without problems. 4. Conclusion Our algorithms are well suited for generating random variates from the (generalized) inverse Gaussian distribution by numerical inversion. In all our experiments the maximal observed u-error was always within acceptable deviation from the requested u-resolution. After a moderate setup, our algorithms are faster than the standard method for generating exponential variates and much faster than numerical root finding algorithms such as the bisection method that has been suggested by Lai (2009). The speed-up factors we observed are between 50 and 220 for sample sizes 105 or larger. Our algorithms are available as C library UNU.RAN and as R package Runuran which can both be downloaded from our web site.
J. Leydold, W. Hörmann / Computational Statistics and Data Analysis 55 (2011) 213–217
217
Acknowledgement The authors are grateful to Carsten Botts for his useful and valuable comments. Appendix. Supplementary data Supplementary material related to this article can be found online at doi:10.1016/j.csda.2010.07.011. References Abramowitz, M., Stegun, I.A. (Eds.), 1972. Handbook of Mathematical Functions, 9th ed.. Dover, New York, ISBN: 0486-612-724. Derflinger, G., Hörmann, W., Leydold, J., 2010. Random variate generation by numerical inversion when only the density is known. ACM Trans. Model. Comput. Simul 20 (4), #18. Galassi, M., et al. GNU Scientific Library. Version 1.13. 2009. http://www.gnu.org/software/gsl/. Hörmann, W., Leydold, J., 2003. Continuous random variate generation by fast numerical inversion. ACM Trans. Model. Comput. Simul. 13 (4), 347–362. Imai, J., Tan, K.S., 2009. an accelerating quasi-Monte Carlo method for option pricing under the generalized hyperbolic Lévy process. SIAM J. Sci. Comput. 31 (3), 2282–2302. doi:10.1137/080727713. Jørgensen, B., 1982. Statistical Properties of the Generalized Inverse Gaussian Distribution. In: Lecture Notes in Statistics, vol. 9. Springer-Verlag, New York, Berlin. Lai, Y., 2009. Generating inverse Gaussian random variates by approximation. Comput. Statist. Data Anal. 53 (10), 3553–3559. doi:10.1016/j.csda.2009.03.007. Leydold, J., Hörmann, W., 2010a. UNU.RAN—a library for non-uniform universal random variate generation. Version 1.6. Department of Statistics and Mathematics, WU Wien, A-1090 Wien, Austria. http://statmath.wu-wien.ac.at/unuran/. Leydold, J., Hörmann, W., 2010b. Runuran—R interface to the UNU.RAN random variate generators. Version 0.13. Department of Statistics and Mathematics, WU Wien, A-1090 Wien, Austria. http://cran.r-project.org/. Neumaier, A., 2001. Introduction to Numerical Analysis. Cambridge University Press, Cambridge. Piessens, R., de Doncker-Kapenga, E., Überhuber, C.W., Kahaner, D.K., 1983. QUADPACK. A Subroutine Package for Automatic Integration. In: Springer Series in Computational Mathematics, vol. 1. Springer-Verlag, Berlin. R Development Core Team. 2009. R: a language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. http://www.R-project.org. ISBN: 3-900051-07-0.