{"id":1243,"date":"2024-06-17T03:00:00","date_gmt":"2024-06-17T01:00:00","guid":{"rendered":"https:\/\/lorentzen.ch\/?p=1243"},"modified":"2025-03-02T18:45:00","modified_gmt":"2025-03-02T17:45:00","slug":"a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution","status":"publish","type":"post","link":"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/","title":{"rendered":"A Tweedie Trilogy \u2014 Part III: From Wrights Generalized Bessel Function to Tweedie&#8217;s Compound Poisson Distribution"},"content":{"rendered":"\n<p><strong>TLDR:<\/strong> The <a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/release.1.7.0.html#scipy-special-improvements\">scipy 1.7.0 release<\/a> introduced Wright\u2019s generalized Bessel function in the Python ecosystem. It is an important ingredient for the density and log-likelihood of <a href=\"https:\/\/en.wikipedia.org\/wiki\/Tweedie_distribution\">Tweedie probabilty distributions<\/a>. In this last part of the trilogy I&#8217;d like to point out why it was important to have this function and share the endeavor of implementing this inconspicuous but highly intractable special function. The fun part is exploiting a free parameter in an integral representation, which can be optimized by curve fitting to the minimal arc length.<\/p>\n\n\n\n<p>This trilogy celebrates the 40th birthday of <a href=\"https:\/\/en.wikipedia.org\/wiki\/Tweedie_distribution\">Tweedie distributions<\/a> in 2024 and highlights some of their very special properties.<\/p>\n\n\n\n<p>See <a href=\"http:\/\/a-tweedie-trilogy-\uff0d-part-i-frequency-and-aggregration-invariance\">part i<\/a> and <a href=\"http:\/\/a-tweedie-trilogy-\uff0d-part-ii-offsets\">part ii<\/a>.<\/p>\n\n\n<div class=\"wp-block-ub-table-of-contents-block ub_table-of-contents\" id=\"ub_table-of-contents-432e28c1-db5d-4329-ae40-7d0f18b066d5\" data-linktodivider=\"false\" data-showtext=\"show\" data-hidetext=\"hide\" data-scrolltype=\"auto\" data-enablesmoothscroll=\"false\" data-initiallyhideonmobile=\"false\" data-initiallyshow=\"true\"><div class=\"ub_table-of-contents-header-container\" style=\"\">\n\t\t\t<div class=\"ub_table-of-contents-header\" style=\"text-align: left; \">\n\t\t\t\t<div class=\"ub_table-of-contents-title\">Table of Contents<\/div>\n\t\t\t\t\n\t\t\t<\/div>\n\t\t<\/div><div class=\"ub_table-of-contents-extra-container\" style=\"\">\n\t\t\t<div class=\"ub_table-of-contents-container ub_table-of-contents-1-column \">\n\t\t\t\t<ul style=\"\"><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#0-tweedie-distributions\" style=\"\">Tweedie Distributions<\/a><\/li><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#1-compound-poisson-gamma\" style=\"\">Compound Poisson Gamma<\/a><\/li><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#2-wrights-generalized-bessel-function\" style=\"\">Wright&#8217;s Generalized Bessel Function<\/a><\/li><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#3-the-integral-representation\" style=\"\">The Integral Representation<\/a><\/li><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#4-arc-length-minimization\" style=\"\">Arc Length Minimization<\/a><\/li><li style=\"\"><a href=\"https:\/\/lorentzen.ch\/index.php\/2024\/06\/17\/a-tweedie-trilogy-part-iii-from-wrights-generalized-bessel-function-to-tweedies-compound-poisson-distribution\/#5-conclusion\" style=\"\">Conclusion<\/a><\/li><\/ul>\n\t\t\t<\/div>\n\t\t<\/div><\/div>\n\n\n<h2 class=\"wp-block-heading\" id=\"0-tweedie-distributions\">Tweedie Distributions<\/h2>\n\n\n\n<p>As pointed out in part I and part II, the family of Tweedie distributions is a very special one with outstanding properties. They are central for estimating expectations with GLMs.  The probability distributions have mainly positive (non-negative) support and are skewed, e.g. Poisson, Gamma, Inverse Gaussian and compound Poisson-Gamma.<\/p>\n\n\n\n<p>As members of the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Exponential_dispersion_model\">exponential dispersion family<\/a>, a slight extension of the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Exponential_family\">exponential family<\/a>, the probability density can be written as<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\nf(y; \\theta, \\phi) &amp;= c(y, \\phi) \\exp\\left(\\frac{y\\theta - \\kappa(\\theta)}{\\phi}\\right)\n\\\\\n\\kappa(\\theta) &amp;= \\kappa_p(\\theta) = \\frac{1}{2-p}((1-p)\\theta)^{\\frac{2-p}{1-p}}\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>It is often more instructive to parametrise the distribution with <code><span class=\"katex-eq\" data-katex-display=\"false\">p<\/span><\/code>, <code><span class=\"katex-eq\" data-katex-display=\"false\">\\mu<\/span><\/code> and <code><code><span class=\"katex-eq\" data-katex-display=\"false\">\\phi<\/span><\/code><\/code>, using<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\n\\theta &amp;= \\begin{cases}\n\\frac{\\mu^{1-p}}{1-p}\\,,\\quad p\\neq 1\\\\\n\\log(\\mu)\\,,\\quad p=1\n\\end{cases}\n\\\\\n\\kappa(\\theta) &amp;= \\begin{cases}\n\\frac{\\mu^{2-p}}{2-p}\\,,\\quad p\\neq 2\\\\\n\\log(\\mu)\\,,\\quad p=2\n\\end{cases}\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>and write<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\nY &amp;\\sim \\mathrm{Tw}_p(\\mu, \\phi)\n\\end{align*}<\/pre><\/div>\n\n\n\n<figure class=\"wp-block-image size-full\"><img loading=\"lazy\" decoding=\"async\" width=\"550\" height=\"464\" src=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image.png\" alt=\"\" class=\"wp-image-1603\" srcset=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image.png 550w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-300x253.png 300w\" sizes=\"auto, (max-width: 550px) 100vw, 550px\" \/><figcaption class=\"wp-element-caption\">Probability density of several Tweedie distributions.<\/figcaption><\/figure>\n\n\n\n<h2 class=\"wp-block-heading\" id=\"1-compound-poisson-gamma\">Compound Poisson Gamma<\/h2>\n\n\n\n<p>A very special domain for the power parameter is between Poisson and Gamma: <code><span class=\"katex-eq\" data-katex-display=\"false\">1&lt;p&lt;2<\/span><\/code>. This range results in the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Compound_Poisson_distribution\">Compound Poisson distribution<\/a> which is suitable if you have a random count process and if each count itself has a random amount. A well know example is insurance claims. Typically, there is a random number of insurance claims, and each and every claim has a random amount of claim costs.<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\nN &amp;\\sim \\mathrm{Poisson}(\\lambda)\\\\\nX_i &amp;\\sim \\mathrm{Gamma}(a, b)\\\\\nY &amp;= \\sum_{i=0}^N X_i \\sim \\mathrm{CompPois}(\\lambda, a, b)\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>For Poisson count we have <code><span class=\"katex-eq\" data-katex-display=\"false\">\\operatorname{E}[N]=\\lambda<\/span><\/code> and <code><span class=\"katex-eq\" data-katex-display=\"false\">\\operatorname{Var}[N]=\\lambda=\\operatorname{E}[N]<\/span><\/code>, for the Gamma amount <code><code><span class=\"katex-eq\" data-katex-display=\"false\">\\operatorname{E}[X]=\\frac{a}{b}<\/span><\/code><\/code> and <code><span class=\"katex-eq\" data-katex-display=\"false\">\\operatorname{Var}[X]=\\frac{a}{b^2}=\\frac{1}{a}\\operatorname{E}[X]^2<\/span><\/code>. For the compound Poisson-Gamma variable, we obtain<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\n\\operatorname{E}[Y] &amp;= \\operatorname{E}[N] \\operatorname{E}[X] = \\lambda\\frac{a}{b}=\\mu\\\\\n\\operatorname{Var}[Y] &amp;=  \\operatorname{Var}[N] \\operatorname{E}[X]^2 +  \\operatorname{E}[N] \\operatorname{Var}[X] =  \\phi \\mu^p\\\\\np &amp;= \\frac{a + 2}{a+1} \\in (1, 2)\\\\\n\\phi &amp;= \\frac{(\\lambda a)^{1-p}}{(p-1)b^{2-p}}\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>What&#8217;s so special here is that there is a point mass at zero, i.e., <code><span class=\"katex-eq\" data-katex-display=\"false\">P(Y=0)=\\exp(-\\frac{\\mu^{2-p}}{\\phi(2-p)}) &gt; 0<\/span><\/code>. Hence, it is a suitable distribution for non-negative quantities with some exact zeros.<\/p>\n\n\n\n<figure class=\"wp-block-image size-full\"><img loading=\"lazy\" decoding=\"async\" width=\"694\" height=\"485\" src=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-1.png\" alt=\"\" class=\"wp-image-1605\" srcset=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-1.png 694w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-1-300x210.png 300w\" sizes=\"auto, (max-width: 694px) 100vw, 694px\" \/><figcaption class=\"wp-element-caption\">Probability density for compound Poisson Gamma, point masses at zero are marked as points.<\/figcaption><\/figure>\n\n\n<div class=\"wp-block-ub-content-toggle wp-block-ub-content-toggle-block\" id=\"ub-content-toggle-block-14991d46-d393-4503-82bc-448601adc5ef\" data-mobilecollapse=\"true\" data-desktopcollapse=\"true\" data-preventcollapse=\"false\" data-showonlyone=\"false\">\n<div class=\"wp-block-ub-content-toggle-accordion\" style=\"border-color: #f1f1f1; \" id=\"ub-content-toggle-panel-block-37a9a858-f583-44ad-991e-9ca030ced861\">\n\t\t\t<div class=\"wp-block-ub-content-toggle-accordion-title-wrap\" style=\"background-color: #f1f1f1;\" aria-controls=\"ub-content-toggle-panel-0-14991d46-d393-4503-82bc-448601adc5ef\" tabindex=\"0\">\n\t\t\t<p class=\"wp-block-ub-content-toggle-accordion-title ub-content-toggle-title-14991d46-d393-4503-82bc-448601adc5ef\" style=\"color: #000000; \">Code<\/p>\n\t\t\t<div class=\"wp-block-ub-content-toggle-accordion-toggle-wrap right\" style=\"color: #000000;\"><span class=\"wp-block-ub-content-toggle-accordion-state-indicator wp-block-ub-chevron-down\"><\/span><\/div>\n\t\t<\/div>\n\t\t\t<div role=\"region\" aria-expanded=\"false\" class=\"wp-block-ub-content-toggle-accordion-content-wrap ub-hide\" id=\"ub-content-toggle-panel-0-14991d46-d393-4503-82bc-448601adc5ef\">\n\n<p><\/p>\n\n\n\n<div class=\"wp-block-codemirror-blocks-code-block code-block\"><pre class=\"CodeMirror\" data-setting='{\"showPanel\":true,\"languageLabel\":\"language\",\"fullScreenButton\":true,\"copyButton\":true,\"mode\":\"python\",\"mime\":\"text\/x-python\",\"theme\":\"material\",\"lineNumbers\":false,\"styleActiveLine\":false,\"lineWrapping\":false,\"readOnly\":true,\"fileName\":\"Python\",\"language\":\"Python\",\"maxHeight\":\"400px\",\"modeName\":\"python\"}'>import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.special import wright_bessel\n\n\ndef cpg_pmf(mu, phi, p):\n    \"\"\"Compound Poisson Gamma point mass at zero.\"\"\"\n    return np.exp(-np.power(mu, 2 - p) \/ (phi * (2 - p)))\n\ndef cpg_pdf(x, mu, phi, p):\n    \"\"\"Compound Poisson Gamma pdf.\"\"\"\n    if not (1 &lt; p &lt; 2):\n        raise ValueError(\"1 &lt; p &lt; 2 required\")\n    theta = np.power(mu, 1 - p) \/ (1 - p)\n    kappa = np.power(mu, 2 - p) \/ (2 - p)\n    alpha = (2 - p) \/ (1 - p)\n    t = ((p - 1) * phi \/ x)**alpha\n    t \/= (2 - p) * phi\n    a = 1 \/ x * wright_bessel(-alpha, 0, t)\n    return a * np.exp((x * theta - kappa) \/ phi)\n\nfig, axes = plt.subplots(ncols=2, figsize=[6.4 * 1.25, 4.8])\nx = np.linspace(1e-9, 10, 200)\nmu = 2\nfor p in [1.2, 1.5, 1.8]:\n    for i, phi in enumerate([0.5, 2]):\n        axes[i].plot(x, cpg_pdf(x=x, mu=mu, phi=phi, p=p), label=f\"{p=}\")\n        axes[i].scatter(0, cpg_pmf(mu=mu, phi=phi, p=p))\n        axes[i].set_ylim(0, 0.5)\n        axes[i].set_title(f\"{phi=}\")\n        if i &gt; 0:\n            axes[i].legend()\n        else:\n            axes[i].set_ylabel(\"pdf(x)\")\n        axes[i].set_xlabel(\"x\")\nfig.suptitle(\"Tweedie Distributions mu=2\")<\/pre><\/div>\n\n<\/div>\n\t\t<\/div>\n<\/div>\n\n\n<p>The rest of this post is about how to compute the density for this parameter range. The easy part is <code><span class=\"katex-eq\" data-katex-display=\"false\">\\exp\\left(\\frac{y\\theta - \\kappa(\\theta)}{\\phi}\\right)<\/span><\/code> which can be directly implemented. The real obstacle is the term <code><span class=\"katex-eq\" data-katex-display=\"false\">c(y, \\phi)<\/span><\/code> which is given by<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\nc(y, \\phi) &amp;= \\frac{\\Phi(-\\alpha, 0, t)}{y}\n\\\\\n\\alpha &amp;= \\frac{2 - p}{1 - p}\n\\\\\nt &amp;= \\frac{\\left(\\frac{(p - 1)\\phi}{y}\\right)^{\\alpha}}{(2-p)\\phi}\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>This depends on <strong>Wright&#8217;s (generalized Bessel) function<\/strong> <code><span class=\"katex-eq\" data-katex-display=\"false\">\\Phi(a, b, z)<\/span><\/code> as introduced in a <a href=\"https:\/\/doi.org\/10.1112\/jlms\/s1-8.1.71\">1933 paper by E. Wright<\/a>.<\/p>\n\n\n\n<h2 class=\"wp-block-heading\" id=\"2-wrights-generalized-bessel-function\">Wright&#8217;s Generalized Bessel Function<\/h2>\n\n\n\n<p>According to <a href=\"https:\/\/dlmf.nist.gov\/10.46#E1\">DLMF 10.46<\/a>, the function is defined as<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{equation*}\n\\Phi(a, b, z) = \\sum_{k=0}^{\\infty} \\frac{z^k}{k!\\Gamma(ak+b)}, \\quad a &gt; -1, b \\in R, z \\in C\n\\end{equation*}<\/pre><\/div>\n\n\n\n<p>which converges everywhere because it is an <a href=\"https:\/\/en.wikipedia.org\/wiki\/Entire_function\">entire function<\/a>. We will focus on the positive real axis <code><span class=\"katex-eq\" data-katex-display=\"false\">z=x\\geq 0<\/span><\/code> and the range <code><span class=\"katex-eq\" data-katex-display=\"false\">a\\geq 0, b\\geq 0<\/span><\/code> (note that <code><span class=\"katex-eq\" data-katex-display=\"false\">a=-\\alpha \\in (0,\\infty)<\/span><\/code> for <code><span class=\"katex-eq\" data-katex-display=\"false\">1&lt;p&lt;2<\/span><\/code>). For the compound Poisson-Gamma, we even have <code><span class=\"katex-eq\" data-katex-display=\"false\">b=0<\/span><\/code>.<\/p>\n\n\n\n<p>Implementation of such a function as done in <a href=\"https:\/\/docs.scipy.org\/doc\/scipy\/reference\/generated\/scipy.special.wright_bessel.html#scipy-special-wright-bessel\">scipy.stats.wright_bessel<\/a>, even for the restricted parameter range, poses tremendous challenges. The first one is that it has three parameters which is quite a lot. Then the series representation above, for instance, can always be used, but depending on the parameters, it will require a huge amount of terms, particularly for large <code><span class=\"katex-eq\" data-katex-display=\"false\">x<\/span><\/code>. As each term involves the <a href=\"https:\/\/en.wikipedia.org\/wiki\/Gamma_function\">Gamma function<\/a>, this becomes expensive very fast. One ends up using different representations and strategies for different parameter regions:<\/p>\n\n\n\n<ul class=\"wp-block-list\">\n<li>Small <code><span class=\"katex-eq\" data-katex-display=\"false\">x<\/span><\/code>: Taylor series according to definition<\/li>\n\n\n\n<li>Small <code><span class=\"katex-eq\" data-katex-display=\"false\">a<\/span><\/code>: Taylor series in <code><span class=\"katex-eq\" data-katex-display=\"false\">a=0<\/span><\/code><\/li>\n\n\n\n<li>Large <code><span class=\"katex-eq\" data-katex-display=\"false\">x<\/span><\/code>: Asymptotic series due to Wright (1935)<\/li>\n\n\n\n<li>Large <code><span class=\"katex-eq\" data-katex-display=\"false\">a<\/span><\/code>: Taylor series according to definition for a few terms around the approximate maximum term <code><span class=\"katex-eq\" data-katex-display=\"false\">k_{max}<\/span><\/code> due to Dunn &amp; Smyth (2005)<\/li>\n\n\n\n<li>General: Integral represantation due to Luchko (2008)<\/li>\n<\/ul>\n\n\n\n<p>Dunn &amp; Smyth investigated several evaluation strategies for the simpler Tweedie density which amounts to Wright&#8217;s functions with <code><span class=\"katex-eq\" data-katex-display=\"false\">b=0<\/span><\/code>, see Dunn &amp; Smyth (2005). Luchko (2008) lists most of the above strategies for the full Wright&#8217;s function.<\/p>\n\n\n\n<p>Note that Dunn &amp; Smyth (2008) provide another strategy to evaluate the Tweedie distribution function by means of the inverse Fourier transform. This does not involve Wright&#8217;s function, but also encounters complicated numerical integration of oscillatory functions. <\/p>\n\n\n\n<h2 class=\"wp-block-heading\" id=\"3-the-integral-representation\">The Integral Representation<\/h2>\n\n\n\n<p>This brings us deep into complex analysis: We start with Hankel&#8217;s contour integral representation of the reciprocal Gamma function.<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{equation*}\n\\frac{1}{\\Gamma(z)} = \\frac{1}{2\\pi i} \\int_{Ha^-} \\zeta^{-z} e^\\zeta \\; d\\zeta\n\\end{equation*}<\/pre><\/div>\n\n\n\n<p>with the Hankel path <code><span class=\"katex-eq\" data-katex-display=\"false\">Ha^-<\/span><\/code> from negative infinity (A) just below the real axis, counter-clockwise with radius <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon&gt;0<\/span><\/code> around the origin and just above the real axis back to minus infinity (D).<\/p>\n\n\n\n<figure class=\"wp-block-image size-full is-resized\"><img loading=\"lazy\" decoding=\"async\" width=\"416\" height=\"256\" src=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/Hankel_contour.drawio.png\" alt=\"\" class=\"wp-image-1639\" style=\"width:650px;height:auto\" srcset=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/Hankel_contour.drawio.png 416w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/Hankel_contour.drawio-300x185.png 300w\" sizes=\"auto, (max-width: 416px) 100vw, 416px\" \/><figcaption class=\"wp-element-caption\">Hankel contour Ha<sup>&#8211;<\/sup> in the complex plane.<\/figcaption><\/figure>\n\n\n\n<p>In principle, one is free to choose <em>any<\/em> such path with the same start (A) and end point (D) as long as one does not cross the negative real axis. One usually lets the AB and CD be infinitesimal close to the negative real line. <strong>Very importantly, the radius <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon&gt;0<\/span><\/code> is a free parameter! That is real magic\ud83e\ude84<\/strong><\/p>\n\n\n\n<p>By interchanging sum and integral and using the series of the exponential, Wright&#8217;s function becomes<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\n\\Phi(a, b, z) &amp;= \\sum_{k=0}^{\\infty} \\frac{z^k}{k!} \\frac{1}{2\\pi i} \\int_{Ha^-} \\zeta^{-(ak+b)} e^\\zeta \\; d\\zeta\n\\\\\n&amp;= \\frac{1}{2\\pi i} \\int_{Ha^-} \\zeta^{-b} e^{\\zeta + z\\zeta^{-a}} \\; d\\zeta\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>Now, one needs to do the tedious work and split the integral into the 3 path sections AB, BC, CD. Putting AB and CD together gives an integral over K, the circle BC gives an integral over P:<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{align*}\n\\Phi(a, b, x) &amp;= \\frac{1}{\\pi} \\int_{\\epsilon}^\\infty K(a, b, x, r) \\; dr\n\\\\\n &amp;+ \\frac{\\epsilon^{1-b}}{\\pi} \\int_0^\\pi P(\\epsilon, a, b, x, \\varphi) \\; d\\varphi\n\\\\\nK(a, b, x, r) &amp;= r^{-b}\\exp(-r + x  r^{-a} \\cos(\\pi a)) \n\\\\\n&amp;\\quad \\sin(x \\cdot r^{-a} \\sin(\\pi a) + \\pi b)\n\\\\\nP(\\epsilon, a, b, x, \\varphi) &amp;= \\exp(\\epsilon \\cos(\\varphi) + x  \\epsilon^{-a}\\cos(a \\varphi))\n\\\\\n&amp;\\quad \\cos(\\epsilon \\sin(\\varphi) - x \\cdot \\epsilon^{-a} \\sin(a \\varphi) + (1-b) \\varphi)\n\\end{align*}<\/pre><\/div>\n\n\n\n<p>What remains is to carry out the numerical integration, also known as quadrature. While this is an interesting topic in its own, let&#8217;s move to the magic part.<\/p>\n\n\n\n<h2 class=\"wp-block-heading\" id=\"4-arc-length-minimization\">Arc Length Minimization<\/h2>\n\n\n\n<p>If you have come so far and say, wow, puh, uff, crazy, \ud83e\udd2f\ud83d\ude31 Just keep on a little bit because here comes the <strong>real fun part<\/strong>\ud83e\udd1e<\/p>\n\n\n\n<p>It turns out that most of the time, the integral over P is the most difficult. The worst behaviour an integrand can have is widely oscillatory. Here is one of my favorite examples:<\/p>\n\n\n\n<figure class=\"wp-block-image size-full\"><img loading=\"lazy\" decoding=\"async\" width=\"623\" height=\"472\" src=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-2.png\" alt=\"\" class=\"wp-image-1642\" srcset=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-2.png 623w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-2-300x227.png 300w\" sizes=\"auto, (max-width: 623px) 100vw, 623px\" \/><figcaption class=\"wp-element-caption\">Integrands for a=5, b=1, x=100 and two choices of epsilon.<\/figcaption><\/figure>\n\n\n\n<p>With the naive choice of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon=1<\/span><\/code>, both integrands (blue) are\u2014well\u2014crazy. There is basically no chance the most sophisticated quadrature rule will work. And then look at the other choice of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon\\approx 4<\/span><\/code>. Both curves seem well behaved (for P, we would need a closer look).<\/p>\n\n\n\n<p>So the idea is to find a good choice of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code> to make P well behaved. Well behaved here means most boring, if possible a straight line. What makes a straight line unique? In flat space, it is the shortest path between two points. Therefore, well behaved integrands have minimal arc length. That is what we want to minimize.<\/p>\n\n\n\n<p>The <a href=\"https:\/\/en.wikipedia.org\/wiki\/Arc_length\">arc length<\/a> <code><span class=\"katex-eq\" data-katex-display=\"false\">S<\/span><\/code> from <code><span class=\"katex-eq\" data-katex-display=\"false\">x=a<\/span><\/code> to <code><span class=\"katex-eq\" data-katex-display=\"false\">x=b<\/span><\/code> of a 1-dimensional function <code><span class=\"katex-eq\" data-katex-display=\"false\">f<\/span><\/code> is given by<\/p>\n\n\n\n<div class=\"wp-block-katex-display-block katex-eq\" data-katex-display=\"true\"><pre>\\begin{equation*}\nS = \\int_a^b \\sqrt{1 + f^\\prime(x)^2} \\; dx\n\\end{equation*}<\/pre><\/div>\n\n\n\n<p>Instead of <code><span class=\"katex-eq\" data-katex-display=\"false\">f=P<\/span><\/code>, we only take the oscillatory part of P and approximate the arc length as <code><span class=\"katex-eq\" data-katex-display=\"false\">f(\\varphi)=f(\\varphi) = \\epsilon \\sin(\\varphi) - x \\epsilon^{-\\rho} \\sin(\\rho \\varphi) + (1-\\beta) \\varphi<\/span><\/code>. For a single parameter point <code><span class=\"katex-eq\" data-katex-display=\"false\">a, b, z<\/span><\/code> this looks like<\/p>\n\n\n\n<figure class=\"wp-block-image size-large\"><img loading=\"lazy\" decoding=\"async\" width=\"1024\" height=\"517\" src=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-3-1024x517.png\" alt=\"\" class=\"wp-image-1648\" srcset=\"https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-3-1024x517.png 1024w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-3-300x151.png 300w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-3-768x387.png 768w, https:\/\/lorentzen.ch\/wp-content\/uploads\/2024\/05\/image-3.png 1211w\" sizes=\"auto, (max-width: 1024px) 100vw, 1024px\" \/><figcaption class=\"wp-element-caption\">Arc length and integrand P for different epsilon, given a=0.1, b=5, x=100.<\/figcaption><\/figure>\n\n\n\n<p>Note the logarithmic y-scale for the right plot of P. The optimal <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon=10<\/span><\/code> is plotted in red and behaves clearly better than smaller values of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code>.<\/p>\n\n\n\n<p>What remains to be done for an actual implementation is<\/p>\n\n\n\n<ul class=\"wp-block-list\">\n<li>Calculate minimal <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code> for a large grid of values <code><span class=\"katex-eq\" data-katex-display=\"false\">a, b, x<\/span><\/code>.<\/li>\n\n\n\n<li>Choose a function with some parameters.<\/li>\n\n\n\n<li>Curve fitting (so again optimisation): Fit this function to the minimal <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code> of the grid via minimising least squares.<\/li>\n\n\n\n<li>Implement some quadrature rules and use this choice of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code> in the hope that it intra- and extrapolates well.<\/li>\n<\/ul>\n\n\n\n<p>This strategy turns out to work well in practice and is implemented in scipy. As the parameter space of 3 variables is huge, the integral representation breaks down in certain areas, e.g. huge values of <code><span class=\"katex-eq\" data-katex-display=\"false\">\\epsilon<\/span><\/code> where the integrands just overflow numerically (in 64-bit floating point precision). But we have our other evaluation strategies for that.<\/p>\n\n\n\n<h2 class=\"wp-block-heading\" id=\"5-conclusion\">Conclusion<\/h2>\n\n\n\n<p>An extensive notebook for Wright&#8217;s function, with all implementation strategies can be found <a href=\"https:\/\/github.com\/lorentzenchr\/notebooks\/blob\/master\/blogposts\/2024-06-17%20wright_bessel_function.ipynb\">here<\/a>.<\/p>\n\n\n\n<p>After an adventurous journey, we arrived at one implementation strategy of Wright&#8217;s generalised Bessel function, namely the integral representation. The path went deep into complex analysis and contour integration, then further to the arc length of a function and finally curve fitting via optimisation. I am really astonished how connected all those different areas of mathematics can be.<\/p>\n\n\n\n<p>Wright&#8217;s function is the missing piece to compute full likelihoods and probability functions of the Tweedie distribution family and is now available in the Python ecosystem via scipy.<\/p>\n\n\n\n<p>We are at the very end of this Tweedie trilogy. I hope it has been entertaining and it has become clear why Tweedie deserves to be celebrated.<\/p>\n\n\n\n<p><strong>Further references:<\/strong><\/p>\n\n\n\n<ul class=\"wp-block-list\">\n<li>Delong, \u0141., Lindholm, M. &amp; W\u00fcthrich, M.V. &#8220;Making Tweedie\u2019s compound Poisson model more accessible&#8221;. <em>Eur. Actuar. J.<\/em> <strong>11<\/strong>, 185\u2013226 (2021). <a href=\"https:\/\/doi.org\/10.1007\/s13385-021-00264-3\">https:\/\/doi.org\/10.1007\/s13385-021-00264-3<\/a><\/li>\n\n\n\n<li>Wright E.M. 1933. &#8220;On the coefficients of power series having essential singularities&#8221;. J. London Math. Soc. 8: 71\u201379. <a href=\"https:\/\/doi.org\/10.1112\/jlms\/s1-8.1.71\">https:\/\/doi.org\/10.1112\/jlms\/s1-8.1.71<\/a><\/li>\n\n\n\n<li>Wright, E. M. 1935, &#8220;The asymptotic expansion of the generalized Bessel&#8221;, function. Proc. London Math. Soc. (2) 38, pp. 257\u2013270. <a href=\"https:\/\/doi.org\/10.1112\/plms\/s2-38.1.257\" target=\"_blank\" rel=\"noreferrer noopener\">https:\/\/doi.org\/10.1112\/plms\/s2-38.1.257<\/a><\/li>\n\n\n\n<li>Dunn, P.K., Smyth, G.K. &#8220;Series evaluation of Tweedie exponential dispersion model densities&#8221;. Stat Comput 15, 267\u2013280 (2005). <a href=\"https:\/\/doi.org\/10.1007\/s11222-005-4070-y\" target=\"_blank\" rel=\"noreferrer noopener\">https:\/\/doi.org\/10.1007\/s11222-005-4070-y<\/a><\/li>\n\n\n\n<li>Dunn, P.K., Smyth, G.K. &#8220;Evaluation of Tweedie exponential dispersion model densities by Fourier inversion&#8221;. Stat Comput 18, 73\u201386 (2008). <a href=\"https:\/\/doi.org\/10.1007\/s11222-007-9039-6\">https:\/\/doi.org\/10.1007\/s11222-007-9039-6<\/a><\/li>\n\n\n\n<li>Luchko, Y. F. (2008), &#8220;Algorithms for Evaluation of the Wright Function for the Real Arguments&#8217; Values&#8221;, Fractional Calculus and Applied Analysis 11(1). <a href=\"https:\/\/eudml.org\/doc\/11309\"> https:\/\/eudml.org\/doc\/11309<\/a><br>Note a slight misprint in the integrand P.<\/li>\n<\/ul>\n","protected":false},"excerpt":{"rendered":"<p>This trilogy celebrates the 40th birthday of Tweedie distributions in 2024 and highlights some of their very special properties.<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[16,17,9],"tags":[6,22],"class_list":["post-1243","post","type-post","status-publish","format-standard","hentry","category-machine-learning","category-programming","category-statistics","tag-python","tag-tweedie-trilogy"],"featured_image_src":null,"author_info":{"display_name":"Christian Lorentzen","author_link":"https:\/\/lorentzen.ch\/index.php\/author\/christian\/"},"_links":{"self":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/1243","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/comments?post=1243"}],"version-history":[{"count":87,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/1243\/revisions"}],"predecessor-version":[{"id":1848,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/posts\/1243\/revisions\/1848"}],"wp:attachment":[{"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/media?parent=1243"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/categories?post=1243"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/lorentzen.ch\/index.php\/wp-json\/wp\/v2\/tags?post=1243"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}