From: Taylor R Campbell Date: Fri, 16 Aug 2019 02:54:44 +0000 (+0000) Subject: Uniform code and style for plots. X-Git-Tag: mit-scheme-pucked-10.1.20~11^2~83 X-Git-Url: https://birchwood-abbey.net/git?a=commitdiff_plain;h=6c8b7f8dd39f7b4a5cbba9232269217e37b14f62;p=mit-scheme.git Uniform code and style for plots. Tweak line widths a little bit to roughly match cmmi10 (Computer Modern Math Italic 10pt) rule widths for axes, and a little thicker for the plots themselves, for the printed manual. --- diff --git a/doc/ref-manual/fig/cn-expm1.eps b/doc/ref-manual/fig/cn-expm1.eps index e64e5e775..8fafaf121 100644 --- a/doc/ref-manual/fig/cn-expm1.eps +++ b/doc/ref-manual/fig/cn-expm1.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,21 +191,22 @@ clip 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)) @@ -179,6 +231,49 @@ clip 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 (plogistic(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} @@ -204,90 +299,19 @@ clip % 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) @@ -295,6 +319,88 @@ clip 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 => --- @@ -328,7 +434,7 @@ clip /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 @@ -339,7 +445,11 @@ clip 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 @@ -350,7 +460,7 @@ clip /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 @@ -392,8 +502,8 @@ clip % 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 @@ -402,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -442,7 +552,7 @@ grestore % 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) diff --git a/doc/ref-manual/fig/cn-log1mexp.eps b/doc/ref-manual/fig/cn-log1mexp.eps index 2fdd6bce1..cd7b822e9 100644 --- a/doc/ref-manual/fig/cn-log1mexp.eps +++ b/doc/ref-manual/fig/cn-log1mexp.eps @@ -74,6 +74,40 @@ clip } 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 (plogistic(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] @@ -153,14 +187,15 @@ clip 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 @@ -192,6 +227,42 @@ clip 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 (plogistic(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) @@ -359,7 +430,7 @@ clip /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 @@ -370,7 +441,11 @@ clip 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 @@ -381,7 +456,7 @@ clip /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 @@ -423,8 +498,8 @@ clip % 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 @@ -433,7 +508,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -465,7 +540,7 @@ grestore gsave 1 0 0 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap 5 dict begin /Tf { diff --git a/doc/ref-manual/fig/cn-log1p.eps b/doc/ref-manual/fig/cn-log1p.eps index cf21a1260..d0294a0cc 100644 --- a/doc/ref-manual/fig/cn-log1p.eps +++ b/doc/ref-manual/fig/cn-log1p.eps @@ -14,8 +14,7 @@ clip % 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 @@ -39,6 +38,10 @@ clip % 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 @@ -54,7 +57,55 @@ clip /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 (plogistic(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. @@ -136,21 +187,22 @@ clip 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)) @@ -175,6 +227,49 @@ clip 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 (plogistic(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} @@ -200,90 +295,19 @@ clip % 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) @@ -291,6 +315,88 @@ clip 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 => --- @@ -324,7 +430,7 @@ clip /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 @@ -335,7 +441,11 @@ clip 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 @@ -346,7 +456,7 @@ clip /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 @@ -354,7 +464,7 @@ clip 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 } { @@ -388,8 +498,8 @@ clip % 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 @@ -398,7 +508,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -437,7 +547,7 @@ grestore gsave 1 0 0 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline -0.9 axmax { dup dup % Tx Tx Tx diff --git a/doc/ref-manual/fig/cn-log1pexp.eps b/doc/ref-manual/fig/cn-log1pexp.eps index 06b4e7b05..346c8b364 100644 --- a/doc/ref-manual/fig/cn-log1pexp.eps +++ b/doc/ref-manual/fig/cn-log1pexp.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,21 +191,22 @@ clip 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)) @@ -179,6 +231,49 @@ clip 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 (plogistic(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} @@ -204,90 +299,19 @@ clip % 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) @@ -295,6 +319,88 @@ clip 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 => --- @@ -328,7 +434,7 @@ clip /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 @@ -339,7 +445,11 @@ clip 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 @@ -350,7 +460,7 @@ clip /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 @@ -392,8 +502,8 @@ clip % 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 @@ -406,7 +516,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont % encoding => encoding { @@ -447,7 +557,7 @@ grestore 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)) diff --git a/doc/ref-manual/fig/cn-logistic.eps b/doc/ref-manual/fig/cn-logistic.eps index 0f5a30ef2..c6f128881 100644 --- a/doc/ref-manual/fig/cn-logistic.eps +++ b/doc/ref-manual/fig/cn-logistic.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,21 +191,22 @@ clip 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)) @@ -179,6 +231,49 @@ clip 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 (plogistic(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} @@ -204,90 +299,19 @@ clip % 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) @@ -295,6 +319,88 @@ clip 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 => --- @@ -328,7 +434,7 @@ clip /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 @@ -339,7 +445,11 @@ clip 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 @@ -350,7 +460,7 @@ clip /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 @@ -392,8 +502,8 @@ clip % 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 @@ -402,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -444,7 +554,7 @@ grestore 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) diff --git a/doc/ref-manual/fig/cn-logistichalf.eps b/doc/ref-manual/fig/cn-logistichalf.eps index 9a8fc476c..4a2ba6cdd 100644 --- a/doc/ref-manual/fig/cn-logistichalf.eps +++ b/doc/ref-manual/fig/cn-logistichalf.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,14 +191,15 @@ clip 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 @@ -179,6 +231,49 @@ clip 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 (plogistic(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} @@ -204,90 +299,19 @@ clip % 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) @@ -295,6 +319,88 @@ clip 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 => --- @@ -328,7 +434,7 @@ clip /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 @@ -339,7 +445,11 @@ clip 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 @@ -350,7 +460,7 @@ clip /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 @@ -392,8 +502,8 @@ clip % 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 @@ -402,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -434,7 +544,7 @@ grestore 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) diff --git a/doc/ref-manual/fig/cn-logit.eps b/doc/ref-manual/fig/cn-logit.eps index 3061a2f64..31f6c515e 100644 --- a/doc/ref-manual/fig/cn-logit.eps +++ b/doc/ref-manual/fig/cn-logit.eps @@ -95,6 +95,19 @@ clip } 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] @@ -174,14 +187,15 @@ clip 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 @@ -213,6 +227,13 @@ clip 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 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) @@ -415,13 +449,32 @@ clip /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 @@ -432,7 +485,11 @@ clip 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 @@ -443,7 +500,7 @@ clip /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 @@ -493,8 +550,8 @@ clip % 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 @@ -524,7 +581,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -559,7 +616,7 @@ grestore 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 diff --git a/doc/ref-manual/fig/cn-logitexp.eps b/doc/ref-manual/fig/cn-logitexp.eps index a60a52535..660586af3 100644 --- a/doc/ref-manual/fig/cn-logitexp.eps +++ b/doc/ref-manual/fig/cn-logitexp.eps @@ -187,7 +187,7 @@ clip 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 @@ -195,7 +195,7 @@ clip 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 @@ -227,6 +227,13 @@ clip 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 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 @@ -434,7 +441,11 @@ clip 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 @@ -487,8 +498,8 @@ clip % 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 @@ -501,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -540,7 +551,7 @@ grestore 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) diff --git a/doc/ref-manual/fig/cn-logithalf.eps b/doc/ref-manual/fig/cn-logithalf.eps index 45c6cfa7d..c8f64488f 100644 --- a/doc/ref-manual/fig/cn-logithalf.eps +++ b/doc/ref-manual/fig/cn-logithalf.eps @@ -74,6 +74,27 @@ clip } 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 (plogistic(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 { @@ -166,14 +187,15 @@ clip 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 @@ -205,16 +227,39 @@ clip 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 (plogistic(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 @@ -396,7 +441,11 @@ clip 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 @@ -449,8 +498,8 @@ clip % 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 @@ -463,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -498,7 +547,7 @@ grestore gsave 1 0 0 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline -.499 .499 { dup dup % Tx Tx Tx diff --git a/doc/ref-manual/fig/cn-loglogistic.eps b/doc/ref-manual/fig/cn-loglogistic.eps index 6fb83712a..2ccccffd8 100644 --- a/doc/ref-manual/fig/cn-loglogistic.eps +++ b/doc/ref-manual/fig/cn-loglogistic.eps @@ -74,6 +74,27 @@ clip } 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 (plogistic(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 { @@ -166,7 +187,7 @@ clip 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 @@ -174,7 +195,7 @@ clip 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 @@ -206,16 +227,39 @@ clip 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 (plogistic(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 @@ -397,7 +441,11 @@ clip 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 @@ -450,8 +498,8 @@ clip % 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 @@ -464,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -504,7 +552,7 @@ grestore gsave 1 0 0 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin axmax { dup Tneg Texp % Tx T(e^-x) diff --git a/doc/ref-manual/fig/expm1.eps b/doc/ref-manual/fig/expm1.eps index 5659358de..485cde35e 100644 --- a/doc/ref-manual/fig/expm1.eps +++ b/doc/ref-manual/fig/expm1.eps @@ -187,7 +187,7 @@ clip 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 @@ -195,7 +195,7 @@ clip 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 @@ -227,6 +227,13 @@ clip 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 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 @@ -491,8 +498,8 @@ clip % 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 @@ -501,7 +508,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -540,7 +547,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin axmax {Texpm1} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/log1mexp.eps b/doc/ref-manual/fig/log1mexp.eps index 78fdc5979..69c351fa8 100644 --- a/doc/ref-manual/fig/log1mexp.eps +++ b/doc/ref-manual/fig/log1mexp.eps @@ -187,7 +187,7 @@ clip 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 @@ -195,7 +195,7 @@ clip 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 @@ -227,6 +227,13 @@ clip 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 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 @@ -490,7 +497,7 @@ clip % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -522,7 +529,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin -.001 {Texp Tneg Tln1p} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/log1p.eps b/doc/ref-manual/fig/log1p.eps index b958b2976..298aec684 100644 --- a/doc/ref-manual/fig/log1p.eps +++ b/doc/ref-manual/fig/log1p.eps @@ -187,7 +187,7 @@ clip 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 @@ -195,7 +195,7 @@ clip 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 @@ -227,6 +227,13 @@ clip 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 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 @@ -491,7 +498,7 @@ clip % Draw reference asymptote. gsave 0.5 setgray - 0.25 setlinewidth + 0.4 setlinewidth [2 3] 0 setdash newpath -1 aymin a2bb moveto @@ -501,7 +508,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -538,7 +545,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline -0.99 axmax {Tln1p} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/log1pexp.eps b/doc/ref-manual/fig/log1pexp.eps index 1bde1cb5f..1dd21bc8f 100644 --- a/doc/ref-manual/fig/log1pexp.eps +++ b/doc/ref-manual/fig/log1pexp.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,21 +191,22 @@ clip 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)) @@ -179,6 +231,49 @@ clip 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 (plogistic(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} @@ -204,90 +299,19 @@ clip % 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) @@ -295,6 +319,88 @@ clip 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 => --- @@ -328,7 +434,7 @@ clip /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 @@ -339,7 +445,11 @@ clip 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 @@ -350,7 +460,7 @@ clip /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 @@ -392,8 +502,8 @@ clip % 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 @@ -402,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont % encoding => encoding { @@ -443,7 +553,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth nspline axmin axmax {Texp Tln1p} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/logistic.eps b/doc/ref-manual/fig/logistic.eps index 9ef061249..ba46fa6bd 100644 --- a/doc/ref-manual/fig/logistic.eps +++ b/doc/ref-manual/fig/logistic.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,21 +191,22 @@ clip 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)) @@ -186,6 +238,42 @@ clip 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 (plogistic(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} @@ -211,90 +299,19 @@ clip % 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) @@ -302,6 +319,88 @@ clip 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 => --- @@ -335,7 +434,7 @@ clip /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 @@ -346,7 +445,11 @@ clip 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 @@ -357,7 +460,7 @@ clip /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 @@ -399,8 +502,8 @@ clip % 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 @@ -409,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -449,7 +552,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin axmax {Tlogistic} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/logistichalf.eps b/doc/ref-manual/fig/logistichalf.eps index 47e03114f..0914ca4db 100644 --- a/doc/ref-manual/fig/logistichalf.eps +++ b/doc/ref-manual/fig/logistichalf.eps @@ -18,8 +18,7 @@ clip % 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 @@ -43,6 +42,10 @@ clip % 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 @@ -58,7 +61,55 @@ clip /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 (plogistic(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. @@ -140,14 +191,15 @@ clip 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 @@ -181,9 +233,45 @@ clip /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 (plogistic(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. @@ -211,90 +299,19 @@ clip % 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) @@ -302,6 +319,88 @@ clip 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 => --- @@ -403,8 +502,8 @@ clip % 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 @@ -417,7 +516,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -449,7 +548,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin axmax {Tlogistic 0.5 Tconst Tsub} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/logit.eps b/doc/ref-manual/fig/logit.eps index cdcd2f953..c60b51008 100644 --- a/doc/ref-manual/fig/logit.eps +++ b/doc/ref-manual/fig/logit.eps @@ -95,6 +95,19 @@ clip } 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] @@ -174,14 +187,15 @@ clip 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 @@ -213,6 +227,13 @@ clip 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 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) @@ -421,7 +455,7 @@ clip /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 @@ -432,7 +466,11 @@ clip 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 @@ -443,7 +481,7 @@ clip /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 @@ -485,8 +523,8 @@ clip % 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 @@ -495,7 +533,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -529,7 +567,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline 0.00001 0.999 {Tlogit} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/logitexp.eps b/doc/ref-manual/fig/logitexp.eps index feadeab86..0a789c1ca 100644 --- a/doc/ref-manual/fig/logitexp.eps +++ b/doc/ref-manual/fig/logitexp.eps @@ -187,7 +187,7 @@ clip 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 @@ -195,7 +195,7 @@ clip 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 @@ -227,6 +227,13 @@ clip 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 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 @@ -434,7 +441,11 @@ clip 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 @@ -487,8 +498,8 @@ clip % 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 @@ -497,7 +508,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -536,7 +547,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin -.001 {Texp Tlogit} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/logithalf.eps b/doc/ref-manual/fig/logithalf.eps index 48778cca8..08f86fb77 100644 --- a/doc/ref-manual/fig/logithalf.eps +++ b/doc/ref-manual/fig/logithalf.eps @@ -74,6 +74,27 @@ clip } 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 (plogistic(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 { @@ -166,14 +187,15 @@ clip 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 @@ -205,16 +227,39 @@ clip 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 (plogistic(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 @@ -396,7 +441,11 @@ clip 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 @@ -449,8 +498,8 @@ clip % 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 @@ -463,7 +512,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -498,7 +547,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline -.499 .499 {Tlogithalf} cubicsplinterpolate grestore diff --git a/doc/ref-manual/fig/loglogistic.eps b/doc/ref-manual/fig/loglogistic.eps index 4e0b611a4..90cf46a94 100644 --- a/doc/ref-manual/fig/loglogistic.eps +++ b/doc/ref-manual/fig/loglogistic.eps @@ -74,6 +74,27 @@ clip } 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 (plogistic(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 { @@ -166,7 +187,7 @@ clip 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 @@ -174,7 +195,7 @@ clip 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 @@ -206,16 +227,39 @@ clip 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 (plogistic(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 @@ -228,9 +272,9 @@ clip /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. @@ -404,7 +448,11 @@ clip 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 @@ -457,8 +505,8 @@ clip % 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 @@ -467,7 +515,7 @@ grestore % Draw axes. gsave - 0.125 setlinewidth + 0.4 setlinewidth /Times-Roman-Numeric /Times-Roman findfont { dup 45 /minus put } reencodefont % replace hyphen by minus @@ -506,7 +554,7 @@ grestore gsave 0 0 1 setrgbcolor - 0.5 setlinewidth + 0.6 setlinewidth 1 setlinecap nspline axmin axmax {Tlogistic Tln} cubicsplinterpolate grestore