% Interpolation parameters.
/nspline 20 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 1 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
% Plot another one: expm1 condition number.
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax { % Tx => T(f(x))
dup Texpm1 % Tx T(e^x-1)
} 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
+
% Automatic differentiation.
/T { 2 array astore } def % x dx => [x dx]
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
/Tlerp % Tt fxmin fxmax => Tx
{
1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
5 dict begin
/Tf {
% Interpolation parameters.
/nspline 20 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
grestore
3 -1 roll add 2 div % tx y ys llx urx (ury+lly)/2
0 exch sub % tx y ys llx urx -(ury+lly)/2
- 3 1 roll sub % tx y ys -(ury+lly)/2 llx-urx
+ 3 1 roll exch sub % tx y ys -(ury+lly)/2 urx-llx
4 index 0 lt {
ticku 2 mul add
} {
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
-1 aymin a2bb moveto
-1 aymax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline -0.9 axmax {
dup dup % Tx Tx Tx
% Interpolation parameters.
/nspline 20 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont % encoding => encoding
{
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
nspline axmin axmax { % Tx => T(f(x))
dup Texp % Tx T(e^x)
dup Tln1p % Tx T(e^x) T(log(1+e^x))
% Interpolation parameters.
/nspline 8 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax { % Tx => T(f(x))
dup Tneg Texp % Tx T(e^-x)
% Interpolation parameters.
/nspline 8 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin 1 a2bb moveto
axmax 1 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax { % Tx => T(f(x))
dup Tneg Texp % Tx T(e^-x)
} 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
+
% Automatic differentiation.
/T { 2 array astore } def % x dx => [x dx]
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
/Tlerp % Tt fxmin fxmax => Tx
{
1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
/ticklabel % v xa ya dxbb dybb => ---
{
4 2 roll a2bb % v dxbb dybb xbb ybb
- newpath moveto rlineto show
+ newpath moveto rmoveto %show
+ gsave
+ dup % v v
+ true charpath pathbbox % v llx lly urx ury
+ 3 -1 roll % v llx urx lly ury
+ 2 copy add 2 div % v llx urx lly ury avgy
+ 3 1 roll % v llx urx avgy lly ury
+ exch sub % v llx urx avgy dy
+ 4 2 roll % v avgy dy llx urx
+ 2 copy add 2 div % v avgy dy llx urx avgx
+ 3 1 roll % v avgy dy avgx llx urx
+ exch sub % v avgy dy avgx dx
+ exch 4 1 roll % v avgx avgy dy dx
+ min % v avgx avgy majaxis
+ 0 360 arc % v
+ closepath
+ 1 setgray
+ fill
+ grestore
+ show
} def
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
0.5 aymin a2bb moveto
0.5 aymax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline 0.00001 0.499 {f} cubicsplinterpolate
nspline 0.501 0.999 {f} cubicsplinterpolate
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin 1 a2bb moveto
axmax 1 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin -0.694 {
dup Texpm1 Tneg % Tx T(1-e^x)
} 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 {
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
-0.5 aymin a2bb moveto
-0.5 aymax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline -.499 .499 {
dup dup % Tx Tx Tx
} 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 {
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin 1 a2bb moveto
axmax 1 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
1 0 0 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax {
dup Tneg Texp % Tx T(e^-x)
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin -1 a2bb moveto
axmax -1 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax {Texpm1} cubicsplinterpolate
grestore
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin -.001 {Texp Tneg Tln1p} cubicsplinterpolate
grestore
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
+ 0.4 setlinewidth
[2 3] 0 setdash
newpath
-1 aymin a2bb moveto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline -0.99 axmax {Tln1p} cubicsplinterpolate
grestore
% Interpolation parameters.
/nspline 20 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont % encoding => encoding
{
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
nspline axmin axmax {Texp Tln1p} cubicsplinterpolate
grestore
% Interpolation parameters.
/nspline 20 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
- Tqv % x dx
+/Texp % Tx => T(e^x)
+{
+ Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
- 1 index mul % e^x e^x*dx
- T % T(e^x)
+ 1 index mul % e^x e^x*dx
+ T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
- Tqv % x dx
- exch dup expm1 % dx x e^x-1
- 3 1 roll % e^x-1 dx x
- 2.718281828 exch exp % e^x-1 dx e^x
- mul % e^x-1 e^x*dx
- T % T(e^x-1)
+ Tqv % x dx
+ exch dup expm1 % dx x e^x-1
+ 3 1 roll % e^x-1 dx x
+ 2.718281828 exch exp % e^x-1 dx e^x
+ mul % e^x-1 e^x*dx
+ T % T(e^x-1)
} def
/Tsin % Tx => T(sin(x))
1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
} def
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
+} def
+
% Interpolation points.
/linpoint { div } def % i n => s_{n,i}
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin 1 a2bb moveto
axmax 1 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax {Tlogistic} cubicsplinterpolate
grestore
% Interpolation parameters.
/nspline 8 def
-/linearsplinterpoint { linpoint } def % i n => s_{i,n}
-/cubicsplinterpoint { chebypoint } def % i n => s_{i,n}
+/splinterpoint { chebypoint } def % i n => s_{i,n}
% Derived axis parameters.
/awidth axmax axmin sub def
% 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
+
/ln1p % x => log(1+x)
{
dup 1.0 add % x 1+x
/expm1
{
- 2.718281828 exch exp 1 sub
+ 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
% Automatic differentiation.
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
/Tlogistic % Tx => T(1/(1+e^{-x}))
{
- 0.0 Tconst exch Tsub Texp % T(e^{-x})
- 1.0 Tconst Tadd % T(1+e^{-x})
- 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
+/Tlerp % Tt fxmin fxmax => Tx
+{
+ 1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
+ 3 -1 roll Tmul % fxmin T((fxmax-fxmin)*t)
+ exch Tconst Tadd % T(fxmin+(fxmax-fxmin)*t)
} def
% Interpolation points.
% Spline interpolation.
-/linearsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0 1 n {/i exch def
- i n linearsplinterpoint Tconst
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb Tq exch Tq exch
- i 0 eq { moveto } { lineto } ifelse
- } for
- stroke
- end
-} def
-
-/cubicsplinterpolate % n fxmin fxmax Tf => ---
-{
- 6 dict begin
- /Tf exch def
- /fxmax exch def
- /fxmin exch def
- /n exch def
- /fxwidth fxmax fxmin sub def
- newpath
- 0.0 0 n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0
- 1 index Tq % Tx0 Ty0 x0
- 1 index Tq % Tx0 Ty0 x0 y0
- moveto % Tx0 Ty0
- 0 1 n 1 sub {/i exch def % Tx0 Ty0
- % Construct a cubic curve segment:
- % c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
- % Note that
- % c(0) = p0, c'(0) = 3 p1 - 3 p0,
- % c(1) = p1, c'(1) = 3 p3 - 3 p2,
- % so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
- % p1 = p0 + c'(0)/3,
- % p2 = p3 - c'(1)/3.
- %
- % Compute Tx3 and Ty3.
- 1.0 i n cubicsplinter
- fxwidth Tconst Tmul fxmin Tconst Tadd
- dup Tf Ta2bb % Tx0 Ty0 Tx3 Ty3
- 4 2 roll % Tx3 Ty3 Tx0 Ty0
- % Compute x1 = x1 + dx1/3.
- exch Tqv % Tx3 Ty3 Ty0 x0 dx0
- 3 div add % Tx3 Ty3 Ty0 x1
- % Compute y1 = y0 + dy0/3.
- exch Tqv % Tx3 Ty3 x1 y0 dy0
- 3 div add % Tx3 Ty3 x1 y1
- % Compute x2 = x3 - dx3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
- 3 div sub % Tx3 Ty3 x1 y1 x2
- % Compute y2 = y3 - dy3/3.
- 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
- 3 div sub % Tx3 Ty3 x1 y1 x2 y2
- % Draw the curve.
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3
- 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
- curveto % Tx3 Ty3
- } for
- pop pop % ---
- stroke
- end
-} def
-
% Compute the tangent vector of a piecewise linear interpolation
% between two consecutive spline interpolation nodes:
%
% c_i(t) := xmin + (xmax - xmin)*(x_i + (x_{i+1} - x_i)*t)
% c_i(1) = c_{i-1}(0)
%
-/cubicsplinter % t i n => T(c_i(t))
+/splinter % t i n => T(c_i(t))
{
2 copy % t i n i n
- cubicsplinterpoint % t i n x0
+ splinterpoint % t i n x0
3 1 roll % t x0 i n
exch 1 add exch % t x0 i+1 n
- cubicsplinterpoint % t x0 x1
+ splinterpoint % t x0 x1
1 index sub Tconst % t x0 T(x1-x0)
3 -1 roll Tvar % x0 T(x1-x0) Tt
Tmul % x0 T((x1-x0)*t)
Tadd % T(x0+(x1-x0)*t)
} def
+% Evaluate the tangent vector f(fxmin + (fxmax - fxmin)*c_i(t)).
+/splinterpoval % fxmin fxmax Tf t i n => x0 y0
+{
+ splinter % fxmin fxmax Tf Tc
+ 4 2 roll % Tf Tc fxmin fxmax
+ Tlerp % Tf Tx0a
+ dup 3 -1 roll % Tx0a Tx0a Tf
+ exec % Tx0a Ty0a
+ Ta2bb % Tx0 Ty0
+} def
+
+% Compute control points of a cubic spline matching the starting and
+% ending tangent vectors. Pops the starting vectors, keeps the ending.
+%
+% c(t) = (1-t)^3 p0 + 3 t (1-t)^2 p1 + 3 t^2 (1-t) p2 + t^3 p3,
+%
+% Note that
+%
+% c(0) = p0, c'(0) = 3 p1 - 3 p0,
+% c(1) = p1, c'(1) = 3 p3 - 3 p2,
+%
+% so p0 = c(0) = (x0, y0), p3 = c(1) = (x1, y1), and
+%
+% p1 = p0 + c'(0)/3,
+% p2 = p3 - c'(1)/3.
+%
+/cubicsplinterpocontrol % Tx0 Ty0 Tx3 Ty3 => Tx3 Ty3 x1 y1 x2 y2
+{
+ 4 2 roll % Tx3 Ty3 Tx0 Ty0
+ % Compute x1 = x1 + dx1/3.
+ exch Tqv % Tx3 Ty3 Ty0 x0 dx0
+ 3 div add % Tx3 Ty3 Ty0 x1
+ % Compute y1 = y0 + dy0/3.
+ exch Tqv % Tx3 Ty3 x1 y0 dy0
+ 3 div add % Tx3 Ty3 x1 y1
+ % Compute x2 = x3 - dx3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x3 dx3
+ 3 div sub % Tx3 Ty3 x1 y1 x2
+ % Compute y2 = y3 - dy3/3.
+ 3 index Tqv % Tx3 Ty3 x1 y1 x2 y3 dy3
+ 3 div sub % Tx3 Ty3 x1 y1 x2 y2
+} def
+
+/cubicsplinterpostart % n fxmin fxmax Tf => Tx0 Ty0 x0 y0
+{
+ 0.0 0 % n fxmin fxmax Tf t i
+ 6 -1 roll % fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0
+ 1 index Tq 1 index Tq % Tx0 Ty0 x0 y0
+} def
+
+/cubicsplinterpostep % Tx0 Ty0 i n fxmin fxmax Tf
+ % => Tx3 Ty3 x1 y1 x2 y2 x3 y3
+{
+ 1.0 % Tx0 Ty0 i n fxmin fxmax Tf t
+ 6 -2 roll % Tx0 Ty0 fxmin fxmax Tf t i n
+ splinterpoval % Tx0 Ty0 Tx3 Ty3
+ cubicsplinterpocontrol % Tx3 Ty3 x1 y1 x2 y2
+ 5 index Tq 5 index Tq % Tx3 Ty3 x1 y1 x2 y2 x3 y3
+} def
+
+/cubicsplinterpostop % Tx3 Ty3 => ---
+{
+ pop pop
+} def
+
+/cubicsplinterpolate % n fxmin fxmax Tf => ---
+{
+ 5 dict begin
+ /Tf exch def
+ /fxmax exch def
+ /fxmin exch def
+ /n exch def
+ /fxwidth fxmax fxmin sub def
+ newpath
+ n fxmin fxmax {Tf} cubicsplinterpostart moveto % Tx0 Ty0
+ 0 1 n 1 sub { n fxmin fxmax {Tf} cubicsplinterpostep curveto } for
+ cubicsplinterpostop
+ stroke
+ end
+} def
+
% Graphics.
/arrowhead % x y angle => ---
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin 0.5 a2bb moveto
axmax 0.5 a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax {Tlogistic 0.5 Tconst Tsub} cubicsplinterpolate
grestore
} 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
+
% Automatic differentiation.
/T { 2 array astore } def % x dx => [x dx]
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+{
+ dup Tq abs 0.5 1.0 3.718281828 div sub le {
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ } {
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ } ifelse
+} def
+
/Tlerp % Tt fxmin fxmax => Tx
{
1 index sub Tconst % Tt fxmin T(fxmax-fxmin)
/xticklabel % x ty => ---
{
exch dup % ty x x
- 3 string cvs % ty x xs
+ 4 string cvs % ty x xs
gsave
newpath 0 0 moveto
dup true charpath
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
/yticklabel % y tx => ---
{
exch dup % tx y y
- 3 string cvs % tx y ys
+ 4 string cvs % tx y ys
gsave
newpath 0 0 moveto
dup true charpath
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
1 aymin a2bb moveto
1 aymax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline 0.00001 0.999 {Tlogit} cubicsplinterpolate
grestore
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
/Tlogit % Tp => T(log(p/(1-p)))
{
dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
} ifelse
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin -.001 {Texp Tlogit} cubicsplinterpolate
grestore
} 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 {
T % T(log(1+x))
} def
-/Texp {
+/Texp % Tx => T(e^x)
+{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
1 index mul % e^x e^x*dx
T % T(e^x)
} def
-/Texpm1
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptotes.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
-0.5 aymin a2bb moveto
-0.5 aymax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline -.499 .499 {Tlogithalf} cubicsplinterpolate
grestore
} 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 {
T % T(log(1+x))
} def
-/Texp % Tx => T(e^x)
+/Texp % Tx => T(e^x)
{
Tqv % x dx
exch 2.718281828 exch exp exch % e^x dx
T % T(e^x)
} def
-/Texpm1 % Tx => T(e^x-1)
+/Texpm1 % Tx => T(e^x-1)
{
Tqv % x dx
exch dup expm1 % dx x e^x-1
T % T(cos(x))
} def
-/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
+/Tlogistic % Tx => T(1/(1+e^{-x}))
+{
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+} def
+
+/Tlogit % Tp => T(log(p/(1-p)))
+{
+ dup Tq -1.0 logistic lt % Tp (p<logistic(-1))
+ 1 index Tq 1.0 logistic gt % Tp (p<logistic(-1)) (p>logistic(1))
+ or {
+ dup 1.0 Tconst exch Tsub % Tp T(1-p)
+ Tdiv Tln % T(log(p/(1-p)))
+ } {
+ dup 2.0 Tconst Tmul % Tp T2p
+ 1.0 Tconst exch Tsub % Tp T(1-2p)
+ exch Tdiv % T((1-2p)/p)
+ Tln1p % T(log1p((1-2p)/p))
+ 0.0 Tconst exch Tsub % T(-log1p((1-2p)/p))
+ } ifelse
+} def
+
+/Tlogithalf % Tp => T(log((1/2+p)/(1/2-p)))
{
dup Tq abs 0.5 1.0 3.718281828 div sub le {
- dup 2.0 Tconst Tmul % Tp T(2p)
- exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
- Tdiv Tln1p % T(log(1+2p/(1/2-p)))
+ dup 2.0 Tconst Tmul % Tp T(2p)
+ exch 0.5 Tconst exch Tsub % T(2p) T(1/2-p)
+ Tdiv Tln1p % T(log(1+2p/(1/2-p)))
} {
- dup 0.5 Tconst Tadd % T(p) T(1/2+p)
- exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
- Tdiv Tln % T(log((1/2+p)/(1/2-p)))
+ dup 0.5 Tconst Tadd % T(p) T(1/2+p)
+ exch 0.5 Tconst exch Tsub % T(1/2+p) T(1/2-p)
+ Tdiv Tln % T(log((1/2+p)/(1/2-p)))
} ifelse
} def
/Tlogistic % Tx => T(1/(1+e^{-x}))
{
- 0.0 Tconst exch Tsub Texp % T(e^{-x})
- 1.0 Tconst Tadd % T(1+e^{-x})
- 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
+ 0.0 Tconst exch Tsub Texp % T(e^{-x})
+ 1.0 Tconst Tadd % T(1+e^{-x})
+ 1.0 Tconst exch Tdiv % T(1/(1+e^{-x}))
} def
% Interpolation points.
2 div % ty x xs ury-lly (urx+llx)/2
0 exch sub % ty x xs ury-lly -(urx+llx)/2
exch % ty x xs -(urx+llx)/2 ury-lly
- ticku 2 mul add % ty x xs -w/2 h
+ 4 index 0 lt {
+ ticku 2 mul add
+ } {
+ pop ticku 2 mul
+ } ifelse % ty x xs -w/2 h
5 -1 roll mul % x xs -w/2 h*ty
4 -1 roll % xs -w/2 h*ty x
0 % xs -w/2 h*ty x y
% Draw reference asymptote.
gsave
0.5 setgray
- 0.25 setlinewidth
- [2 3] 0 setdash
+ 0.4 setlinewidth
+ [4 4] 0 setdash
newpath
axmin axmin a2bb moveto
axmax axmax a2bb lineto
% Draw axes.
gsave
- 0.125 setlinewidth
+ 0.4 setlinewidth
/Times-Roman-Numeric
/Times-Roman findfont
{ dup 45 /minus put } reencodefont % replace hyphen by minus
gsave
0 0 1 setrgbcolor
- 0.5 setlinewidth
+ 0.6 setlinewidth
1 setlinecap
nspline axmin axmax {Tlogistic Tln} cubicsplinterpolate
grestore
projects / mit-scheme.git / commitdiff