File: invigamma.dem

package info (click to toggle)
gnuplot 6.0.2%2Bdfsg1-2
  • links: PTS, VCS
  • area: main
  • in suites: sid
  • size: 14,940 kB
  • sloc: ansic: 95,319; cpp: 7,590; makefile: 2,470; javascript: 2,328; sh: 1,531; lisp: 664; perl: 304; pascal: 191; tcl: 88; python: 46
file content (51 lines) | stat: -rw-r--r-- 1,086 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#
# Unit test for inverse incomplete gamma integral
#
# z = invigamma( a, p )
# p = igamma( a, z )
#
# Domain restrictions
# 	0 < p <= 1
#       a > 0

set title "Inverse normalized lower incomplete gamma function"
set xrange [0:1]
set yrange [0:8]
set xlabel "p"
set ylabel "invigamma( a, p )"
set key outside box opaque invert reverse Left
set grid
set tics scale 0
set sample 500

plot invigamma( 0.5, x ) lw 3 title "a = 0.5", \
     invigamma( 1.0, x ) lw 3 title "a = 1.0", \
     invigamma( 1.5, x ) lw 3 title "a = 1.5", \
     invigamma( 2.0, x ) lw 3 title "a = 2.0"

pause -1 "<cr> to continue"
reset 

set title "Residual after igamma(invigamma(a,p))"
set xlabel "p"
set ylabel "residual"
set sample 200

residual(a,p) =  igamma( a, invigamma(a,p) ) - p

set xrange [0:1]

do for [E=10:-2:-1] {
    plot residual( 10**E, x ) title sprintf("a = %g", 10**E)
    pause .5
}

pause -1 "<cr> to continue"

E = -3
set title "Residual after igamma(invigamma(a,p))" .\
          "\nSuffers from numerical overflow/underflow when a < 0.005"
replot

pause -1 "<cr> to continue"
reset