Sign in to follow this  
Followers 0
andybiochem

[SOLVED] Maths problem (lower incomplete Gamma function)

7 posts in this topic

#1 ·  Posted (edited)

Hi,

PROBLEM:

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

Posted Image

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:

Posted Image

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.

#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 by andybiochem

- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites



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 ;-)

#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).


This wonderful site allows debugging and testing regular expressions (many flavors available). An absolute must have in your bookmarks.
Another excellent RegExp tutorial. Don't forget downloading your copy of up-to-date pcretest.exe and pcregrep.exe here
RegExp tutorial: enough to get started
PCRE v8.33 regexp documentation latest available release and currently implemented in AutoIt beta.

SQLitespeed is another feature-rich premier SQLite manager (includes import/export). Well worth a try.
SQLite Expert (freeware Personal Edition or payware Pro version) is a very useful SQLite database manager.
An excellent eBook covering almost every aspect of SQLite3: a must-read for anyone doing serious work.
SQL tutorial (covers "generic" SQL, but most of it applies to SQLite as well)
A work-in-progress SQLite3 tutorial. Don't miss other LxyzTHW pages!
SQLite official website with full documentation (may be newer than the SQLite library that comes standard with AutoIt)

Share this post


Link to post
Share on other sites

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.


- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

#5 ·  Posted (edited)

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.

post-44800-12567693815783_thumb.jpg

You can PM any time you wish.

Edited by jchd

This wonderful site allows debugging and testing regular expressions (many flavors available). An absolute must have in your bookmarks.
Another excellent RegExp tutorial. Don't forget downloading your copy of up-to-date pcretest.exe and pcregrep.exe here
RegExp tutorial: enough to get started
PCRE v8.33 regexp documentation latest available release and currently implemented in AutoIt beta.

SQLitespeed is another feature-rich premier SQLite manager (includes import/export). Well worth a try.
SQLite Expert (freeware Personal Edition or payware Pro version) is a very useful SQLite database manager.
An excellent eBook covering almost every aspect of SQLite3: a must-read for anyone doing serious work.
SQL tutorial (covers "generic" SQL, but most of it applies to SQLite as well)
A work-in-progress SQLite3 tutorial. Don't miss other LxyzTHW pages!
SQLite official website with full documentation (may be newer than the SQLite library that comes standard with AutoIt)

Share this post


Link to post
Share on other sites

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.


This wonderful site allows debugging and testing regular expressions (many flavors available). An absolute must have in your bookmarks.
Another excellent RegExp tutorial. Don't forget downloading your copy of up-to-date pcretest.exe and pcregrep.exe here
RegExp tutorial: enough to get started
PCRE v8.33 regexp documentation latest available release and currently implemented in AutoIt beta.

SQLitespeed is another feature-rich premier SQLite manager (includes import/export). Well worth a try.
SQLite Expert (freeware Personal Edition or payware Pro version) is a very useful SQLite database manager.
An excellent eBook covering almost every aspect of SQLite3: a must-read for anyone doing serious work.
SQL tutorial (covers "generic" SQL, but most of it applies to SQLite as well)
A work-in-progress SQLite3 tutorial. Don't miss other LxyzTHW pages!
SQLite official website with full documentation (may be newer than the SQLite library that comes standard with AutoIt)

Share this post


Link to post
Share on other sites

#7 ·  Posted (edited)

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!!!! :mad:

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. :) :) :P:idea::lol:

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

:D

Edited by andybiochem

- Table UDF - create simple data tables - Line Graph UDF GDI+ - quickly create simple line graphs with x and y axes (uses GDI+ with double buffer) - Line Graph UDF - quickly create simple line graphs with x and y axes (uses AI native graphic control) - Barcode Generator Code 128 B C - Create the 1/0 code for barcodes. - WebCam as BarCode Reader - use your webcam to read barcodes - Stereograms!!! - make your own stereograms in AutoIT - Ziggurat Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Box-Muller Gaussian Distribution RNG - generate random numbers based on normal/gaussian distribution - Elastic Radio Buttons - faux-gravity effects in AutoIT (from javascript)- Morse Code Generator - Generate morse code by tapping your spacebar!

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!


Register a new account

Sign in

Already have an account? Sign in here.


Sign In Now
Sign in to follow this  
Followers 0