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
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
% Math utilities.
/eps 1.0 { dup 2 div dup 1.0 add 1.0 eq { pop exit } if exch pop } loop def
/sqrteps eps sqrt def
/cbrteps eps 1 3 div exp def
/pi 3.14159265358979 def
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
dup 1.0 eq { pop } {
dup ln % x 1+x log(1+x)
3 -1 roll % 1+x log(1+x) x
mul % 1+x x*log(1+x)
exch % x*log(1+x) 1+x
1.0 sub % x*log(1+x) (1+x)-1
div % x*log(1+x)/[(1+x)-1]
} ifelse
} def
/expm1
{
dup abs cbrteps gt {
% e^x - 1
2.718281828 exch exp 1 sub
} { dup abs sqrteps gt {
% x + x^2/2 + x^3/3
dup dup dup dup % x x x x
mul mul 6 div % x x x^3/6
exch dup mul 2 div % x x^3/6 x^2/2
add add % x+x^2/2+x^3/6
} { dup abs eps gt {
% x + x^2/2
dup dup mul 2 div add
} {
% nothing
} ifelse } ifelse } ifelse
} def
/logistic % x => 1/(1+e^{-x})
{
0 exch sub 2.718281828 exch exp 1 add 1 exch div
} def
/logit % p => log(p/(1-p))
{
dup -1.0 logistic lt % p (p<logistic(-1))
1 index 1.0 logistic gt % p (p<logistic(-1)) (p>logistic(1))
or {
dup 1.0 exch sub % p 1-p
div ln % log(p/(1-p))
} {
dup 2.0 mul % p 2p
1.0 exch sub % p 1-2p
exch div % (1-2p)/p
ln1p % log1p((1-2p)/p)
0.0 exch sub % -log1p((1-2p)/p)
} ifelse
} def
/logithalf % p => log((1/2+p)/(1/2-p))
{
dup abs 0.5 1.0 3.718281828 div sub le {
dup 2 mul % p 2p
exch 0.5 exch sub % 2p 1/2-p
div ln1p % log(1+2p/(1/2-p))
} {
dup 0.5 add % p 1/2+p
exch 0.5 exch sub % 1/2+p 1/2-p
div ln % log((1/2+p)/(1/2-p))
} ifelse
} def
/versin % x => versin(x)
{
2.0 div % x/2
sin % sin(x/2)
dup mul % sin^2(x/2)
2.0 mul % 2*sin^2(x/2)
} def
|
