[SOLVED] Maths problem (lower incomplete Gamma function)

[andybiochem]

Hi,

PROBLEM:

I'm trying to write an AI script that will be able to calculate probabilities from the CHI squared cumulative distribution function:

This is calculated from Gamma function and lower incomplete gamma function.

I'm not sure if my code is out, or if it is my interpretation of the maths. My gamma function (Euler's infinate product definition) looks ok - I've checked it against known coefficients. I think the problem may lie with my lower incomplete.

Anyway, this is what I get:

Note the red line where k=1 and chi > ~4. Degrees of freedom greater than 1 seem to work out ok.

Hoping there's some maths guru here on the forum.

expandcollapse popup

#include "BigNum.au3"
#include "GraphGDIPlus UDF.au3"

Opt("GUIOnEventMode", 1)

$GUI = GUICreate("",600,600)
GUISetOnEvent(-3,"_Exit")
GUISetState()

;----- Create Graph area -----
$Graph = _GraphGDIPlus_Create($GUI,40,30,530,520,0xFF000000,0xFFaabbcc)

;----- Set X axis range from -5 to 5 -----
Global $xLow = 0
Global $xHigh = 8
Global $yLow = 0
Global $yHigh = 1
Global $PointsToPlot = 100
_GraphGDIPlus_Set_RangeX($Graph,$xLow,$xHigh,8,1,1)
_GraphGDIPlus_Set_RangeY($Graph,$yLow,$yHigh,10,1,1)

;----- Set Y axis range from -5 to 5 -----
_GraphGDIPlus_Set_GridX($Graph,1,0xFF99aabb,0xFF778899)
_GraphGDIPlus_Set_GridY($Graph,0.1,0xFF99aabb,0xFF778899)

;----- Draw the graph -----
_Draw_Graph()

While 1
    Sleep(100)
WEnd

Func _Draw_Graph()
;----- Put Function To Test Here -----

For $k = 1 to 5
    _GraphGDIPlus_Set_PenColor($Graph,0xFF342356)
    _GraphGDIPlus_Set_PenSize($Graph,1)

    ;----- draw lines -----
    $First = True
    For $x = $xLow to $xHigh Step ($xHigh - $xLow) / $PointsToPlot
        $y = _CHISquaredPDF($k,$x)
        If $First = True Then _GraphGDIPlus_Plot_Start($Graph,$x,$y)
        $First = False
        _GraphGDIPlus_Plot_Line($Graph,$x,$y)
        _GraphGDIPlus_Refresh($Graph)
    Next

Next

EndFunc


Func _CHISquaredPDF($k,$x)

    Return _Lower_IncompleteGammaFunction($k/2,$x/2) / _GammaFunction($k/2)

EndFunc



Func _GammaFunction($iZ)
    $nPCError = 0.00001 ;allowable % error
    $nProduct = 1
    $n = 0
    Do
        $n += 1
        $nProductPre = $nProduct
        $y = ((1 + (1/$n))^$iZ) / (1 + ($iZ / $n))
        $nProduct *= $y
    Until (($nProduct - $nProductPre) / $nProduct) * 100 <= $nPCError And $n > 1
    Return (1/$iZ) * $nProduct
EndFunc



Func _Lower_IncompleteGammaFunction($iS,$iZ)
    $iPCError = 0.00001
    $iSum = 0
    $iSumPre = 0
    $iKFactorial = 1
    $k = 0
    Do
        If $k > 0 Then $iKFactorial = _BigNum_Mul($iKFactorial,$k)
        $iSumPre = $iSum
        $iSum += (((-1)^$k) / ($iKFactorial)) * (($iZ^($iS + $k)) / ($iS + $k))
        $k += 1
    Until Abs((($iSum - $iSumPre) / $iSum) * 100) <= $iPCError Or $k > 1754
    Return $iSum
EndFunc



Func _Exit()
    ;----- close down GDI+ and clear graphic -----
    _GraphGDIPlus_Delete($GUI,$Graph)
    Exit
EndFunc

GraphGDIPlus UDF.au3 is here

BigNum.au3 is here

Edited October 29, 2009 by andybiochem

[PsaltyDS]

Dang, reviewed every episode of Numb3rs on my DVR, and found nothing on it...

[jchd]

I'm not sure if my code is out, or if it is my interpretation of the maths. My gamma function (Euler's infinate product definition) looks ok - I've checked it against known coefficients. I think the problem may lie with my lower incomplete.

Hi,

Your code has several problems, some minor but _one_ bigger.

First, and that's the main reason for the visible issue at k = 1, your Γ($iZ) incorrectly stops _way_ before it should. Euler series converges, but only _very_ slowly. You should stop at $n > 500000 (not 1) or even higher if you whish to get below the accuracy you pretend. This is because you're summing very small terms, getting always closer to the preceding. Thus you reach you precision percentage, but you're far from the _total_ accuracy you're after.

Then, why do you compute _GammaFunction($_k/2) for each increment in ... $x ? That's a big saving.

Also, as a bonus saving, computing (-1)^$k for increasing integer $k can be streamlined to doing an alternate summation.

FYI, what you're computing isn't the PDF (Probability Density Function) but is instead the CDF (Cumulative Distribution Function).

I have double-checked the numerical results against Mathematica v5.1 and there is no significant problem with the attached code. Do not try to push decimals too far anyway. You'd probably easier to precompute Γ(z) for commonly used values of z (with Im(z) = 0) and store them in a table. Also 0.00001 % for fitting curves may be a bit pedantic ;-)

expandcollapse popup

#include "BigNum.au3"
#include "GraphGDIPlus.au3"

Opt("GUIOnEventMode", 1)

$GUI = GUICreate("",600,600)
GUISetOnEvent(-3,"_Exit")
GUISetState()

;----- Create Graph area -----
$Graph = _GraphGDIPlus_Create($GUI,40,30,530,520,0xFF000000,0xFFaabbcc)

;----- Set X axis range from -5 to 5 -----
Global $xLow = 0
Global $xHigh = 8
Global $yLow = 0
Global $yHigh = 1
Global $PointsToPlot = 100
_GraphGDIPlus_Set_RangeX($Graph,$xLow,$xHigh,8,1,1)
_GraphGDIPlus_Set_RangeY($Graph,$yLow,$yHigh,10,1,1)

;----- Set Y axis range from -5 to 5 -----
_GraphGDIPlus_Set_GridX($Graph,1,0xFF99aabb,0xFF778899)
_GraphGDIPlus_Set_GridY($Graph,0.1,0xFF99aabb,0xFF778899)

;----- Draw the graph -----
_Draw_Graph()

While 1
    Sleep(100)
WEnd

Func _Draw_Graph()
;----- Put Function To Test Here -----

For $_k = 1 to 5
    _GraphGDIPlus_Set_PenColor($Graph,0xFF342356)
    _GraphGDIPlus_Set_PenSize($Graph,1)

    ;----- draw lines -----
    $First = True
    $Gamma_k = _GammaFunction($_k/2)
    For $x = $xLow to $xHigh Step ($xHigh - $xLow) / $PointsToPlot
;~         $y = _CHISquaredPDF($_k,$x)      ;; not what we want (name change)
;~         $y = _CHISquaredCDF($_k,$x)      ;; avoid computing Γ($_k) for each point!
        $y = _Lower_IncompleteGammaFunction($_k/2,$x/2) / $Gamma_k
        If $First = True Then _GraphGDIPlus_Plot_Start($Graph,$x,$y)
        $First = False
        _GraphGDIPlus_Plot_Line($Graph,$x,$y)
        _GraphGDIPlus_Refresh($Graph)
    Next

Next

EndFunc


;; This is _not_ the PDF (Probability Density Function)
;
;; It is instead the CDF (Cumulative Distribution Function)
;;
;~ Func _CHISquaredPDF($k,$x)
Func _CHISquaredCDF($k,$x)

    Return _Lower_IncompleteGammaFunction($k/2,$x/2) / _GammaFunction($k/2)

EndFunc



Func _GammaFunction($iZ)
    $nPCError = 0.00001 ;allowable % error
    $nProduct = 1
    $n = 0
    Do
        $n += 1
        $nProductPre = $nProduct
        $y = ((1 + (1/$n))^$iZ) / (1 + ($iZ / $n))
        $nProduct *= $y
    Until (($nProduct - $nProductPre) / $nProduct) * 100 <= $nPCError And $n > 300000
_ConsoleWrite('Γ(' & $iZ & ') = ' & (1/$iZ) * $nProduct & @LF)
    Return (1/$iZ) * $nProduct
EndFunc



Func _Lower_IncompleteGammaFunction($iS,$iZ)
    $iPCError = 0.00001
    $iSum = 0
    $iSumPre = 0
    $iKFactorial = 1
    $SumK = 0
    $Sign = 1
    Do
        If $Sumk > 0 Then $iKFactorial = _BigNum_Mul($iKFactorial,$SumK)
        $iSumPre = $iSum
        $iSum += $Sign / ($iKFactorial) * (($iZ^($iS + $SumK)) / ($iS + $SumK))
        $Sign *= -1
        $SumK += 1
    Until Abs((($iSum - $iSumPre) / $iSum) * 100) <= $iPCError Or $SumK > 1754
;~ _ConsoleWrite('γ(' & $iS & ', ' & $iZ & ') = ' & $iSum & @LF)
    Return $iSum
EndFunc



Func _ConsoleWrite($sString)
    Local $aResult = DllCall("kernel32.dll", "int", "WideCharToMultiByte", "uint", 65001, "dword", 0, "wstr", $sString, "int", -1, _
                                "ptr", 0, "int", 0, "ptr", 0, "ptr", 0)
    If @error Then Return SetError(1, @error, 0)
    Local $tText = DllStructCreate("char[" & $aResult[0] & "]")
    $aResult = DllCall("Kernel32.dll", "int", "WideCharToMultiByte", "uint", 65001, "dword", 0, "wstr", $sString, "int", -1, _
                            "ptr", DllStructGetPtr($tText), "int", $aResult[0], "ptr", 0, "ptr", 0)
    If @error Then Return SetError(2, @error, 0)
    ConsoleWrite(DllStructGetData($tText, 1))
EndFunc


Func _Exit()
    ;----- close down GDI+ and clear graphic -----
    _GraphGDIPlus_Delete($GUI,$Graph)
    Exit
EndFunc

Note: save this as UTF-8 source to have greek symbols show up (sorry, bad habit of mine).

[andybiochem]

Thanks a MILLION jchd. I'm no mathmatician...or rather, a very slooow learning one.

First, and that's the main reason for the visible issue at k = 1, your Γ($iZ) incorrectly stops _way_ before it should. Euler series converges, but only _very_ slowly. You should stop at $n > 500000 (not 1) or even higher if you whish to get below the accuracy you pretend.

I understand. The $n > 1 term was only there to prevent the Do ending on the first loop (which it was doing), not as an attempt at accuracy. I was hoping my PCError would end the loop at a suitably accurate point. Now I know I was way out.

Then, why do you compute _GammaFunction($_k/2) for each increment in ... $x ? That's a big saving.

Ha! I didn't spot that! Thanks again!

FYI, what you're computing isn't the PDF (Probability Density Function) but is instead the CDF (Cumulative Distribution Function).

Oops, my bad. I knew I was computing CDF. I had previously used the same script to compute PDF, and forgot to change the func names.

You'd probably easier to precompute Γ(z) for commonly used values of z (with Im(z) = 0) and store them in a table.

I really wanted to get away from that approach. I converted a c++ math version of gamma,lower incomplete gamma, complemented, inverse, log normal etc etc... but all these functions used look-up coefficients. I want to be able to calculate these on the fly with 'some' accuracy.

Also 0.00001 % for fitting curves may be a bit pedantic ;-)

Heh! This example was just to see if I was getting the right values (which I wasn't). I'll be using these functions in some bigger stats programs I'm developing.

Thanks again jchd, I really appreciate your help.

In truth, it's so rarely that I come across anyone I can discuss this stuff with. I work in a hospital laboratory, so much of my time is spent with practical analysis rather than fundamental maths. Yet, this sort of stuff is very important to what I do: normal distribution, chi squared...and as a result gamma function etc.

As a non-mathmatician, am I attacking this sort of calculus in the right way? Are there any books you'd recommend that cover maths (especially stats) in programming??

Much respect.

[jchd]

Hi,

Glad it could help you.

To tell you the truth, I'm not at all a mathematician either, nor a statistician, indeed nor anything else than a slave working like crazy without getting ever paid ;-)

I nonetheless have kept some interest in math, but more toward number theory and other specialized topics like integral lattices, [ternary] quadratic forms and other unfriendly beasts no normal human cares about. I really can't help you much with statistics and to be fully honest, I sometimes regard them as highly suspiscuous. Of course I'm kind of joking here and I understand that they can be of paramount interest in lab environment.

I have no idea what other simple algorithm you can use for Gamma having better convergence characteristic. Google around and you will probably find useful ideas. Another possibility, if you have to routinely perform such computations, is to make a wrapper for a decently supported open-source math library. This way you gain the speed and robustness of a proven library without loosing the flexibility of AutoIt.

Sorry, I can't recommend any textbook on this kind of subject but again Google is your best bet. For your eyes only, here's a Mma one-liner.

You can PM any time you wish.

Edited October 28, 2009 by jchd

[jchd]

If you insist on having your code free of external libraries, you may want to check some of the algorithms presented on this page as it appears they could well fit your bill without eating too many cycles nor depending on "magic" constants.

I didn't try any but just by browsing the page I found some formulas amazingly simple yet very rapidly converging.

Since this is now getting way off topic, I suggest to shift further discussion to PM or mail (easier), except if other AutoIt users claim interest.

[andybiochem]

Thanks for the link!

I've adjusted my code, and am well on my way to finishing that section of my stats prog.

Again, thanks.

EDIT

Hold on one minute!!!!

I've just noticed I'm missing an ABS() around the PCError calculation in my OP Gamma Function:

Func _GammaFunction($iZ)
    $nPCError = 0.00001 ;allowable % error
    $nProduct = 1
    $n = 0
    Do
        $n += 1
        $nProductPre = $nProduct
        $y = ((1 + (1/$n))^$iZ) / (1 + ($iZ / $n))
        $nProduct *= $y
    Until (($nProduct - $nProductPre) / $nProduct) * 100 <= $nPCError And $n > 1
    Return (1/$iZ) * $nProduct
EndFunc

This...

Until (($nProduct - $nProductPre) / $nProduct) * 100 <= $nPCError And $n > 1

...should read ...

Until Abs((($nProduct - $nProductPre) / $nProduct) * 100) <= $nPCError And $n > 1

Because we're converging around the 'true' value, the PCError calulation sometimes comes out negative, which means it satisfies <= 0.00001

Thus, the converge ends after two iterations when k=1. I think you spotted this in your original post, but I misunderstood.

Adding the absolution solves the problem with my original script. :)

Anyway, I learnt a lot here! ...and I've made the optimisations you suggested!

Edited October 29, 2009 by andybiochem