gammaUpper_reg_inv — Gamma

The inverse of the regularized upper incomplete Gamma integral

View versions (2)

Interface

#include <codecogs/maths/special/gamma/gammaupper_reg_inv.h>

using namespace Maths::Special::Gamma;

double gammaUpper_reg_inv(double a, double y0)

This function returns the inverse of the regularized "upper" incomplete Gamma Integral, thereby finding the root of:

$$y - f_a(x) = 0$$

(1)

where

$$f_a(x) = \frac{1}{\Gamma(a)} \int_x^{\infty} e^{-t} t^{a-1}\ dt$$

(2)

In other words, given y0, this function finds x such that

$$\mathrm{gammaUpper\_reg}(a,x) = y_0$$

(3)

It therefore finds the lower limit of the Gamma integral.

This function starts with an approximate value of

$$x = a t^3$$

(4)

where $t = 1 - d - CDF_{inv}(y) \sqrt{d}$ and $d = \frac{1}{9a}$

The routine will perform up to 10 Newton iterations to find the root of gammaUpper_reg(a,x) - y = 0.

Accuracy:

Tested at random a, y in the intervals indicated. <pre> Relative error: domain(a,y) domain # trials peak rms 0.5,100 0,0.5 100000 1.0e-14 1.7e-15 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14 </pre>

Example 1

#include <stdio.h>
#include <codecogs/maths/special/gamma/gammaupper_reg_inv.h>
#include <codecogs/maths/special/gamma/gamma_upper_reg.h>

int main()
{
  for(double x=1, a=0.01; x<10; x+=0.5,a+=0.1)
  {
    double y = Maths::Special::Gamma::gammaUpper_reg(a, x);
    printf("\n gammaUpper_reg(%lf,%lf)=%lf",a,x,y);
    double x2= Maths::Special::Gamma::gammaUpper_reg_inv(a, y);
    printf("   gammaUpper_reg_inv(%lf, %lf)=%lf", a,y,x2);
  }
  return 0;
}

Output:

gammaUpper_reg(0.010000,1.000000)=0.002216   gammaUpper_reg_inv(0.010000, 0.002216)=1.000000
 gammaUpper_reg(0.110000,1.500000)=0.012627   gammaUpper_reg_inv(0.110000, 0.012627)=1.500000
 gammaUpper_reg(0.210000,2.000000)=0.013812   gammaUpper_reg_inv(0.210000, 0.013812)=2.000000
 gammaUpper_reg(0.310000,2.500000)=0.012406   gammaUpper_reg_inv(0.310000, 0.012406)=2.500000
 gammaUpper_reg(0.410000,3.000000)=0.010395   gammaUpper_reg_inv(0.410000, 0.010395)=3.000000
 gammaUpper_reg(0.510000,3.500000)=0.008433   gammaUpper_reg_inv(0.510000, 0.008433)=3.500000
 gammaUpper_reg(0.610000,4.000000)=0.006726   gammaUpper_reg_inv(0.610000, 0.006726)=4.000000
 gammaUpper_reg(0.710000,4.500000)=0.005311   gammaUpper_reg_inv(0.710000, 0.005311)=4.500000
 gammaUpper_reg(0.810000,5.000000)=0.004168   gammaUpper_reg_inv(0.810000, 0.004168)=5.000000
 gammaUpper_reg(0.910000,5.500000)=0.003259   gammaUpper_reg_inv(0.910000, 0.003259)=5.500000
 gammaUpper_reg(1.010000,6.000000)=0.002542   gammaUpper_reg_inv(1.010000, 0.002542)=6.000000
 gammaUpper_reg(1.110000,6.500000)=0.001979   gammaUpper_reg_inv(1.110000, 0.001979)=6.500000
 gammaUpper_reg(1.210000,7.000000)=0.001540   gammaUpper_reg_inv(1.210000, 0.001540)=7.000000
 gammaUpper_reg(1.310000,7.500000)=0.001197   gammaUpper_reg_inv(1.310000, 0.001197)=7.500000
 gammaUpper_reg(1.410000,8.000000)=0.000930   gammaUpper_reg_inv(1.410000, 0.000930)=8.000000
 gammaUpper_reg(1.510000,8.500000)=0.000723   gammaUpper_reg_inv(1.510000, 0.000723)=8.500000
 gammaUpper_reg(1.610000,9.000000)=0.000561   gammaUpper_reg_inv(1.610000, 0.000561)=9.000000
 gammaUpper_reg(1.710000,9.500000)=0.000436   gammaUpper_reg_inv(1.710000, 0.000436)=9.500000

Parameters

the main argument, which must be greater than zero.

the target value y such that gammaUpper_reg(a,x) = y0. Must also be greater than zero.

Interactive Calculator

Result –