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
52
53
54
55
56
57
58
59
60
61
62
|
c function pchisq.f
c
c source
c Bevington, pages 192-193.
c
c purpose
c evaluate probability for exceeding chi square
c
c usage
c result = pchisq (chisqr, nfree)
c
c description of parameters
c chisqr - comparison value of reduced chi square
c nfree - number of degrees of freedom
c
c subroutines and function subprograms required
c gamma (x)
c calculates gamma function for integers and half-integers
c
c comments
c calculation is approximate for nfree odd and
c chi square greater than 50
c
function pchisq (chisqr,nfree)
double precision z,term,sum
11 if (nfree) 12,12,14
12 pchisq=0.
goto 60
14 free=nfree
z=chisqr*free/2.
neven=2*(nfree/2)
if (nfree-neven) 21,21,41
c
c number of degrees of freedom is even
c
21 imax=nfree/2
term=1.
sum=0.
31 do 34 i=1,imax
fi=i
sum=sum+term
34 term=term*z/fi
35 pchisq=sum*dexp(z)
goto 60
c
c number of degrees of freedom is odd
c
41 if (z-25) 44,44,42
42 z=chisqr*(free-1.)/2.
goto 21
44 pwr=free/2.
term=1.
sum=term/pwr
51 do 56 i=1,1000
fi=i
term=-term*z/fi
sum=sum+term/(pwr+fi)
55 if (dabs(term/sum)-0.00001) 57,57,56
56 continue
57 pchisq=1.-(z**pwr)*sum/gamma(pwr)
60 return
end
|
