gamma Upper reg inv

The inverse of the regularized upper incomplete Gamma integral

Controller: CodeCogs

Contents

Interface
GammaUpper Reg Inv
Page Comments

Dependents

Interface

C++

Excel

GammaUpper Reg Inv

doublegammaUpper_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$

(2)

where

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

(3)

In otherwords, give y, this fucntion finds x such that

gammaUpper_reg(a,x) = y

It therefore find the lower limit of the Gamma integral.

The function computes the inverse of the regularized upper incomplete Gamma integral, by finding the lower limit for the Gamma integral.

This function, start within an approximate value of

$x = a t^3$

(4)

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

The routine will performs 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

a	the man argument, which must be greater than zero.
y0	the result of equivalence to a call from gammaUpper_reg(a,x). Must also be greater than zero.

Authors

Stephen L.Moshier. Copyright 1984, 1987, 1989, 1992, 2000
Documented and Updated by Will Bateman

Source Code

Source code is available when you agree to a GP Licence or buy a Commercial Licence.

Not a member, then Register with CodeCogs. Already a Member, then Login.