Symmetric Stable PDF (McCulloch 95)

/*************  Symmetric Stable PDF (McCulloch 95)  *********************
Subject: Symmetric Stable PDF (McCulloch)

Posted by: J. Huston McCulloch
           Economics Dept.
           Ohio State Univ.
           1945 N. High St.
           Columbus, OH 43210
           (614) 292-0382
           mcculloch.2@osu.edu

[Revised 9/20/95].
These programs were written by J. Huston McCulloch, June 1993.
The code is written and submitted for public, non-commerical use,
with no performance guarantees.  The computation is described in
J. Huston McCulloch, "Numerical Approximation of the Symmetric Stable
Distributions and Densities," Ohio State Univ. Econ Dept, Oct. 1994,
which is available from the author on request.  Any research use
of this code should cite this working paper.

There are 3 similar programs:
symstb1 returns the symmetric stable density s(x; alpha) only.
symstb2 returns s(x; alpha) along with d/dx s(x; alpha).
symstb3 returns the complemented cdf 1-S(x; alpha), along with
  s(x; alpha) and d/dx s(x; alpha).
two procs, named poly and dgamma, follow symstb3, and are used by
  all 3 programs.

Accuracy:
The expected relative density precision is 1.0e-6 for alpha in the
range [.84, 2.00].  The programs have considerably reduced precision for
alpha < .84, although no error message is given for lower alpha
values.  The absolute precision of the density is 6.6e-5
for alpha in the range [..92, 2.00], while that of the distribution is
2.2e-5 in the same range.  See paper for details.

Speed:
Changing alpha induces set-up calculations, so submit as many x values
as you can before changing alpha.  Submit x in long vectors if possible.
Densities takes an average of 33 microseconds per ordinate
on a P-5/100 using GAUSS 3.2.12 if ordinates are submitted to
symstb1 in repeated vectors of length 1000 with a common alpha value.
Symstb2 and symstb3 are respectively slower.

Compatibility:
The programs use several global variables.  Their names have all been
"uglified" by use of a terminal underscore (_) to make conflicts unlikely.

GAUSS code:
************************************************************/
proc symstb1(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns sden = symmetric stable density @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_;
declare cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden;
local x11, x2a1, x2pp, cd, gd;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sden = gden(x);
  retp(sden);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
oldalf_ = alpha;

startx:
xa1 = 1 + a_ * abs(x);
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 -1) /a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
cd = cden (x1);
gd = gden(x2);
j1 = trunc(20 * z) + 1;
rpz = poly(pd_[.,j1], z, 4);
sden = (alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz) .* zp;
retp (sden);
endp;

proc(2) = symstb2(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns {sden,sdenp} @
@ sden = symmetric stable density @
@ sdenp = derivative of same wrt x @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_;
declare cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xx, xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden;
local sign, zpp, x11, x1pp, x2a1, x2pp, cd, cdp, gd, gdp, rppz, cgr, sdenp;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sden = gden(x);
  sdenp = -sden .* x / 2;
  retp(sden, sdenp);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
pdd_ = pd_[2:5,.] .* seqa (1,1,4);
oldalf_ = alpha;

startx:
sign = 2 * (x .>= zeros(rows(x),1)) - 1;
xx = abs(x);
xa1 = 1 + a_ * xx;
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
zpp = -(alpha + 1) * a_ * zp ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x1pp = 2 * x11 .* x1p;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 - 1) / a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
x2pp = 1.5 * x2p ./ xa1a;
cd = cden (x1);
cdp = - 2 * pi * x1 .* cd .* cd;
gd = gden(x2);
gdp = - x2 .* gd / 2;
j1 = trunc(20 * z) + 1;
rpz = poly(pd_[.,j1], z, 4);
rppz = poly(pdd_[.,j1], z, 3);
cgr = alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz;
sden = cgr .* zp;
sdenp = sign .* (cgr .* zpp + zp .* zp .* (alf2_ * (cd.*x1pp + x1p.*x1p.*cdp)
         + alf1_ * (gd .* x2pp + x2p .* x2p .* gdp) -rppz));
retp (sden, sdenp);
endp;

proc(3) = symstb3(x,alpha);
@ June 1993: GAUSS version @
@ If you use this program, please cite @
@  J. Huston McCulloch, "Numerical Approximation of the Symmetric @
@  Stable Distribution and Density," OSU Econ Dept., Oct. 1994.  @
@ x may be a column vector; alpha is a scalar @
@ returns {sdisc,sden,sdenp} @
@ sdisc = complemented symmetric stable distribution @
@ sden = symmetric stable density @
@ sdenp = derivative of density wrt x @
declare s_, sqpi_, a2_, cpxp0_, gpxp0_, cpxpp0_, gpxpp0_, cppp_, gppp_, c5i_;
declare a_, alf2_, alf1_, pd_, p_, pdd_, sq2_;
declare first_ = 1;
declare oldalf_ = -999.99;
local r, q, znot, zn4, zn5;
local i, i1, j, j1;
local pialf, sp0, sppp0, xp0, xpp0, xppp0, spzp1;
local rp0, rpp0, rppp0, rp1, b, c;
local xx, xa1, xa1a, z, x1, x2, zp, x1p, x2p, rpz, sden, pos, sdisc, rz;
local sign, zpp, x11, x1pp, x2a1, x2pp, cd, cdp, gd, gdp, rppz, cgr, sdenp;
fn afun(alpha) = 2^(1/alpha) - 1;
fn cden(x) = 1/(pi*(1+x.*x));
fn gden(x) = exp(maxc(-700*ones(1,rows(x)) | -(x.*x)'/4)) / (2*sqpi_);
fn cfun(x) = .5 - atan(x) / pi;
fn gfun(x) = cdfnc (x/sq2_);
fn zfun(x, alpha, a) = 1 - (1 + a*x) ^ (-alpha);
fn zpfun(x, alpha, a) = alpha * a * (1+a*x) ^ (-alpha -1);
fn xfun(z, alpha, a) = ((1-z)^(-1/alpha) - 1) / a;

if(not first_); goto starta; endif;
@ this s is the gij of the writeup @
s_ =    {1.8514190959D2,   -4.6769332663D2,
         4.8424720302D2,   -1.7639153404D2,
        -3.0236552164D2,    7.6351931975D2,
        -7.8560342101D2,    2.8426313374D2,
         4.4078923600D2,   -1.1181138121D3,
         1.1548311335D3,   -4.1969666223D2,
        -5.2448142165D2,    1.3224487717D3,
        -1.3555648053D3,    4.8834079950D2,
         5.3530435018D2,   -1.3374570340D3,
         1.3660140118D3,   -4.9286099583D2,
        -4.8988957866D2,    1.2091418165D3,
        -1.2285872257D3,    4.4063174114D2,
         3.2905528742D2,   -7.3211767697D2,
         6.8183641829D2,   -2.2824291084D2,
        -2.1495402244D2,    3.9694906604D2,
        -3.3695710692D2,    1.0905855709D2,
         2.1112581866D2,   -2.7921107017D2,
         1.1717966020D2,    3.4394664342D0,
        -2.6486798043D2,    1.1999093707D2,
         2.1044841328D2,   -1.5110881541D2,
         9.4105784123D2,   -1.7221988478D3,
         1.4087544698D3,   -4.2472511892D2,
        -2.1990475933D3,    4.2637720422D3,
        -3.4723981786D3,    1.0174373627D3,
         3.1047490290D3,   -5.4204210990D3,
         4.2221052925D3,   -1.2345971177D3,
        -5.1408260668D3,    1.1090264364D4,
        -1.0270337246D4,    3.4243449595D3,
         1.1215157876D4,   -2.4243529825D4,
         2.1536057267D4,   -6.8490996103D3,
        -1.8120631586D4,    3.1430132257D4,
        -2.4164285641D4,    6.9126862826D3,
         1.7388413126D4,   -2.2108397686D4,
         1.3397999271D4,   -3.1246611987D3,
        -7.2435775303D3,    4.3545399418D3,
         2.3616155949D2,   -7.6571653073D2,
        -8.7376725439D3,    1.5510852129D4,
        -1.3789764138D4,    4.6387417712D3};
s_ = reshape (s_,19, 4);
sqpi_ = sqrt(pi);
sq2_ = sqrt(2);
a2_ = afun(2);
cpxp0_ = 1/pi;
gpxp0_ = 1/(4*sqpi_*a2_);
cpxpp0_ = 2/pi;
gpxpp0_ = gpxp0_ * 1.5;
cppp_ = 4/pi;
gppp_ = gpxpp0_ * 2.5 - 1/(32*sqpi_*a2_^3);
c5i_ = {1, 5, 10, 10, 5, 1};
first_ = 0;

starta:
if alpha > 2.000001;
  format /rdn 20,15;
  print "alpha too big: " alpha;
  end;
  endif;
if alpha >= 2;
  sdisc = gfun(x);
  sden = gden(x);
  sdenp = -sden .* x / 2;
  retp(sdisc, sden, sdenp);
  endif;
if alpha == oldalf_; goto startx; endif;
znot = seqa(1,1,19) /20;
zn4 = (1-znot)^4;
zn5 = zn4 .* (1-znot);
a_ = afun(alpha);
alf2_ = 2 - alpha;
alf1_ = alpha - 1;
pialf = pi * alpha;
sp0 = dgamma(1/alpha) / pialf;
sppp0 = -dgamma(3/alpha) / pialf;
xp0 = 1/(alpha * a_);
xpp0 = xp0 * (1+alpha) / alpha;
xppp0 = xpp0 * (1+2*alpha)/alpha;
spzp1 = (a_^alpha) * dgamma(alpha) * sin(pialf / 2) / pi;
rp0 = -sp0 * xp0 + alf2_ * cpxp0_ + alf1_ * gpxp0_;
rpp0 = -sp0 * xpp0 + alf2_ * cpxpp0_ + alf1_ * gpxpp0_;
rppp0 = -sp0 * xppp0 - sppp0 * xp0^3 + alf2_ * cppp_ + alf1_ * gppp_;
rp1 = -spzp1 + alf2_ / pi;
q = zeros(6,1);
q[2] = rp0;
q[3] = rpp0/2;
q[4] = rppp0/6;
r = alf2_ * (s_ * (alf2_^seqa(1,1,4) - 1));
b = -sumc(q[2:4]) - sumc(r.*zn5);
c = rp1 - q[2:4]' * seqa(1,1,3) - 5 * sumc(r.*zn4);
q[5] = 5*b-c;
q[6] = c-4*b;
p_ = zeros(6,20);
i = 0;
do until i > 5;
i1 = i + 1;
j = 1;
do until j > 19;
j1 = j + 1;
p_[i1, j1] = p_[i1, j] + r[j] * (-znot[j])^(5-i);
j = j1;
endo;
i = i1;
endo;
p_ = q + c5i_ .* p_;
pd_ = p_[2:6,.] .* seqa (1,1,5);
pdd_ = pd_[2:5,.] .* seqa (1,1,4);
oldalf_ = alpha;

startx:
pos = x .>= zeros(rows(x),1);
sign = 2 * pos - 1;
xx = abs(x);
xa1 = 1 + a_ * xx;
xa1a = xa1 ^ (-alpha);
z = 1 - xa1a;
zp = alpha * a_ * xa1a ./ xa1;
zpp = -(alpha + 1) * a_ * zp ./ xa1;
x11 = 1/xa1a;
x1 = x11 - 1;
x1p = x11 .* x11;
x1pp = 2 * x11 .* x1p;
x2a1 = 1/sqrt(xa1a);
x2 = (x2a1 - 1) / a2_;
x2p = x2a1 ./ (2*a2_ * xa1a);
x2pp = 1.5 * x2p ./ xa1a;
cd = cden (x1);
cdp = - 2 * pi * x1 .* cd .* cd;
gd = gden(x2);
gdp = - x2 .* gd / 2;
j1 = trunc(20 * z) + 1;
rz = poly(p_[.,j1], z, 5);
rpz = poly(pd_[.,j1], z, 4);
rppz = poly(pdd_[.,j1], z, 3);
cgr = alf2_ * cd .* x1p + alf1_ * gd .* x2p - rpz;
sden = cgr .* zp;
sdenp = sign .* (cgr .* zpp + zp .* zp .* (alf2_ * (cd.*x1pp + x1p.*x1p.*cdp)
         + alf1_ * (gd .* x2pp + x2p .* x2p .* gdp) -rppz));
sdisc = alf2_*cfun(x1) + alf1_*gfun(x2) + rz;
sdisc = sdisc .* pos + (1-sdisc) .* (1-pos);
retp (sdisc, sden, sdenp);
endp;

proc poly(a, x, k);
local sum, j;
sum = a[k+1,.]';
j = k;
do until j ==0;
sum = sum .* x + a[j,.]';
j = j - 1;
endo;
retp(sum);
endp;

proc dgamma(x);
@ Note -- GAUSS 3.2+ has adopted an improved gamma function similar
to this one, making this proc redundant. @
local t, ser, c, stp;
c = {76.18009172947146, -86.50532032941677, 24.01409824083091,
     -1.231739572450155, .001208650973866179, -.5395239384953d-5};
stp = 2.5066282746310005;
t = x + 5.5;
t = t^(x+.5)/exp(t);
ser = 1.000000000190015 + sumc(c./(x+seqa(1,1,6)));
retp(t * stp * ser / x);
endp;