{"id":43,"date":"2010-05-14T15:52:13","date_gmt":"2010-05-14T13:52:13","guid":{"rendered":"http:\/\/djalil.chafai.net\/blog\/?p=43"},"modified":"2019-11-10T19:09:17","modified_gmt":"2019-11-10T18:09:17","slug":"simulating-heavy-tails","status":"publish","type":"post","link":"https:\/\/djalil.chafai.net\/blog\/2010\/05\/14\/simulating-heavy-tails\/","title":{"rendered":"Simulating heavy tails"},"content":{"rendered":"

\"\"<\/a>How to quickly simulate a positive real random variable $X$ with a polynomial tail, say $\\mathbb{P}(X>t)\\sim (1+t)^{-\\alpha}$ as $t\\to\\infty$, for some prescribed $\\alpha>0$. We may use an already coded function provided by\u00a0 our favorite computational software. Quoting my former PhD advisor, you just have to press the button rand or something similar<\/em>. But how to proceed from scratch using uniform random numbers? The inversion method suggests to set $X=U^{-1\/\\alpha}-1$ where $U$ is a uniform random variable on $[0,1].$ The law of $X$ is Pareto with probability density function $t\\in[0,\\infty)\\mapsto a(1+t)^{-(1+a)}$.\u00a0 It is amusing to remark that the fractional part of $X$ has probability density function<\/p>\n

$$t\\in[0,1]\\mapsto\\alpha\\sum_{n=1}^\\infty\\frac{1}{(n+t)^{1+\\alpha}}=\\alpha\\zeta(\\alpha+1,t+1)$$<\/p>\n

where $\\zeta$ is the Hurwitz Zeta function<\/a>. The integer part of $X$ follows a Zipf like law.<\/p>\n

To simulate an $n$-sample of $X$ with GNU-Octave<\/a> one can use rand(1,n).^(-1\/alpha)-1<\/em>. The problem is that on computers, real numbers are replaced by\u00a0 floats<\/a>: there is a quantification limit.\u00a0 Small values of $\\alpha$ produce a lack of precision for large values, which are frequent for heavy tails...<\/p>\n

Regarding the algorithm, it is not wise to use the inversion method to simulate a law on an infinite set (recall that there is no uniform law of such sets).\u00a0 To simulate $X$, one should prefer an excess packing method as in Marsaglia & Tsang Ziggurat algorithm<\/a>. This alternative approach has the advantage of being faster and a bit more precise (the code is heavier). It is implemented for exponential, Gaussian, and gamma distributions in the\u00a0 rande<\/em>, randg<\/em>, randn<\/em> functions of GNU-Octave 3.2. If you read French, you may take a look at the first chapter of my pedagogical book with Bernard Bercu<\/a>. This approach constitutes an instance of\u00a0 a general meta principle in Computer Science which states that one can speed up algorithms by using static pre-computations, related to the rigid structure of the problem. The fastest algorithm is not necessarily the shortest algorithm in terms of number of code lines.<\/p>\n

If you are interested in the simulation of $\\alpha$-stable distributions, $0<\\alpha\\leq 2$, you may take a look at the Chambers-Mallows-Stuck exact algorithm<\/a> which generalizes the well known Box-Muller algorithm for Gaussians (suffers from the quantification limit).<\/p>\n

For true (unpredictable) randomness, take a look at random.org<\/strong><\/a> and its GNU-R interface<\/strong><\/a> \ud83d\ude09<\/p>\n","protected":false},"excerpt":{"rendered":"

How to quickly simulate a positive real random variable $X$ with a polynomial tail, say $\\mathbb{P}(X>t)\\sim (1+t)^{-\\alpha}$ as $t\\to\\infty$, for some prescribed $\\alpha>0$. We may…<\/p>\n

Continue readingSimulating heavy tails<\/span><\/a><\/div>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"iawp_total_views":57},"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/43"}],"collection":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/comments?post=43"}],"version-history":[{"count":6,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/43\/revisions"}],"predecessor-version":[{"id":11756,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/posts\/43\/revisions\/11756"}],"wp:attachment":[{"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/media?parent=43"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/categories?post=43"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/djalil.chafai.net\/blog\/wp-json\/wp\/v2\/tags?post=43"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}