/******************* Integrate Multivariate Normal ******************

Author: Geert Molenberghs <gmolenb@luc.ac.be>
Date: 2 Sep 96

This code is provided without guarantee for public non-commercial use.
**********************************************************************/

/*
**                   file MVN.ARC
**
** GAUSS-version of the algorithm by M. Shervish to integrate the multivariate
** normal density.
**
** Geert Molenberghs, April 1990
** Biostatistics, Limburgs Universitair Centrum
** Universitaire Campus, Bldg D
** B 3590 Diepenbeek
** Belgium
** gmolenb@luc.ac.be
**
** The author cannot be held responsible for errors in the code
**
** (The headers of the procedures are in English, unfortunately, the comments
** in the middle of the procedures is in Dutch).
**
*/

/*
** CDFNINV
**
** FORMAT   y=3Dcdfninv(x);
**
** INPUT    x is a value between 0 and 1
**
** OUTPUT   the quantile    =20
**
*/

proc cdfninv(x);
/* z zoeken zodanig dat cdfn(z)=3Dx */

local w1,w2,hulp,i;

/* testen op slechte input */

if ((x<=3D-0.00001) or (x>=3D1.00001));
   errorlog "de input ligt buiten [0,1]";
   end;
endif;

if x<=3D0;
   retp(-10);
   end;
elseif x>=3D1;
   retp(10);
   end;
endif;

/* de eindpunten van het zoekinterval fixeren */

if (cdfn(0)<x);
   w2=3D1;
   do until (cdfn(w2)>=3Dx);
      w2=3Dw2+1;
   endo;
   w1=3Dw2-1;
else;
   w1=3D-1;
   do until (cdfn(w1)<x);
      w1=3Dw1-1;
   endo;
   w2=3Dw1+1;
endif;

/* het zoekinterval verfijnen */

i=3D1;

do while ((abs(w1-w2)>10E-10) and (i<=3D50));

   hulp=3D(w1+w2)/2;
   if (cdfn(hulp)<x);
      w1=3Dhulp;
   else;
      w2=3Dhulp;
   endif;
   i=3Di+1;
endo;
if (i=3D=3D50);
   errorlog "cdfninv heeft een te lage nauwkeurigheid (geen stop)";
endif;
retp(w1);
endp;


/*
** PHI - aux. procedure - not for direct use
**
** FORMAT   t=3Dphi(x,y);
**
** INPUT    x is a real value
**          y value realted to the conditional density.
**
** OUTPUT   t calculated funcional value
**
** REMARK   phi calculates univariate normal density
**
*/

proc phi(x,y);

local arg,sq2p,xlow,result;

/*
** XLOW is een waarde die op mogelijke underflow test. Zo het argument
** onder XLOW ligt wordt de waarde 0 teruggegeven.
*/

xlow=3D-87;

/*
** sq2p geeft de logaritme van de wortel uit 2*pi. Zo kan die mee in de
** exponent verhuizen en daalt het aantal berekeningen.
*/

sq2p=3D0.91893853320467274;
result=3D0;
arg=3D-0.5*x*x-sq2p-y;
if (arg>xlow);
   result=3Dexp(arg);
endif;
retp(result);
endp;


/*
** Nu volgt een reeks functies, ge=8Bmplementeerd als procedures, die moeten
** dienen om afgeleide functies te berekenen met het oog op de fouten-
** afschatting bij de diverse Newton-Cotes regels.
*/

proc f2(x,y);
retp(abs((x*x-1)/(y*y)));
endp;

proc f3(x,y);
retp(abs(-x*(x*x-3)/(y*y*y)));
endp;

proc f4(x,y);
retp(abs((3+x*x*(-6+x*x))/(y^4)));
endp;

proc f5(x,y);
retp(abs(x*(-15+x*x*(10-x*x)))/(y^5));
endp;

proc f6(x,y);
retp(abs(-15+x*x*(45-x*x*(-15+x*x)))/(y^6));
endp;


/*
**  MULVAR - key procedure
**
**  FORMAT  {prob,bound}=3Dmulvar(a,b,sigma,eps,n,inf);
**
**  OUTPUT  prob   : calculated integral (1x1)
**          bound  : calculated error margin
**
**  INPUT    a     : the n upper bounds (nx1)
**           b     : the n lower bounds (nx1)
**           sigma : correlationmatrix   (nxn)
**           eps   : pre-specified error margin(1x1)
**           n     : dimension of the integral (1x1)
**           inf   : vector of 'infinities'   (nx1)
**                      1(0): lower boundary (upper boundary) is infinite
**                      2   : the whole interval is finite
**
**  REMARKS
**
** This procedure is the main one of the implementation in GAUS of the=20
** algorithm by M.J. Shervish:
**  Multivariate Normal Probabilities with Error Bound, Royal Statistical
**  Society (1984). Algorithm AS 195.
*/


proc (2) =3D mulvar(a,b,sigma,eps,n,inf);

/* Definitie van de lokale variabelen */

/* constanten */

local coeff,ddo,eo,cons2,cons3,cons4,epsmin,epssim,bcs,bcn;


/* matrices */

/* bladzijde B */
local bound,ept,z,cup;
/* bladzijde C */
local sinv,dett,lrhpm,tel,cnt,simp,prob,sd,rho,nm1,nm2;
/* bladzijde D */
local c,d,t,eplos,fac,ipr,x1,x2,y1,y2,epsi;
/* bladzijde E */
local l,hm,cond,condl,beta,i,j,k,sigc,binc,bh;
/* bladzijde F */
local co,iend,ien;
/* bladzijde G */
local xss,xm,prod,pr2,ni;
/* bladzijde H */
local bint,intvl,ans,bou4,prep,preb;
/* bladzijde I */
local nl,s,num,ok,njs,mjs,inj,numl;
/* bladzijde J */
local ep,r,fe,ind,bl;
/* bladzijde K */
local br,getal,hulp,bou5,del;
/* bladzijde L */
local hulp1,hulp2,hulp3,bou1,bou2,bou3,itype,ij,jk;
/* bladzijde N */
local ksa;
/* bladzijde P */
local fact,xs,ll;
/* bladzijde Q */
local ans;
/* bladzijde R */
local tem,temp,indl,ity,temb,fsa;
/* bladzijde S */
local intl,igo;

/* initialisatie */

let eo[13,1]=3D  -2.85697     -2.3344142    -1.7320508
               -1.3556262   -1            -0.74196378
                0            0.74196378    1
                1.3556262    1.7320508     2.3344142
                2.85697;

let ddo[25,1]=3D -3.75043972  -3.3242574    -2.85697
               -2.36675941  -2.3344142    -1.8891759
               -1.7320508   -1.3556262    -1.15440539
               -1           -0.74196378   -0.61670659
                0            0.61670659    0.74196378
                1            1.15440539    1.3556262
                1.7320508    1.8891759     2.3344142
                2.36675941   2.85697       3.3242574
                3.75043972;

cons2=3D6.879833;
cons3=3D4.517004;
cons4=3D76.371214;
epsmin=3D1E-8;
epssim=3D6E-5;

let coeff[5,3]=3D 0.3111111111111 0.333333333333 0.5
                1.4222222222222 0              0
                0.5333333333333 1.333333333333 0
                1.4222222222222 0              0
                0.3111111111111 0.333333333333 0.5;

sd=3Dzeros(2,1);
cond=3Dzeros(7,1);
condl=3Dzeros(7,1);
beta=3Dzeros(7,7);
binc=3Dzeros(7,7);
xss=3Dzeros(7,7);
xm=3Dzeros(7,1);
prod=3Dzeros(7,1);
pr2=3Dzeros(7,1);
bint=3Dzeros(7,1);
intvl=3Dzeros(7,1);
ans=3Dzeros(7,1);
bou4=3Dzeros(7,1);
prep=3Dzeros(7,1);
preb=3Dzeros(7,1);
s=3Dzeros(7,36);
num=3Dzeros(7,1);
ep=3Dzeros(7,1);
r=3Dzeros(36,7);
fe=3Dzeros(7,1);
ind=3Dzeros(7,1);
bl=3Dzeros(7,7);
bou5=3Dzeros(7,1);
del=3Dzeros(7,7);
itype=3Dzeros(7,1);
beta=3Dzeros(7,7);

let br=3D br1 br2 br3 br4 br5 br6 br7 ;

tel=3D1;
do while (tel<=3D7);
   cnt=3Dvarput(zeros(7,7),br[tel]);
   tel=3Dtel+1;
endo;

let bcs[5,1]=3D1 4 6 4 1;
let bcn[7,1]=3D1 6 15 20 15 6 1;

ksa=3Dzeros(7,1);
xm=3Dzeros(7,1);
fact=3Dzeros(7,1);
ans=3Dzeros(7,1);

/*
** BLADZIJDE B: foutendetectie van input en bepaling Z-waarde
**
** Via CDFNINV(EPT) berekenen we een waarde X zodat CDFN(x)<EPT.
** De waarde Z vervangt verderop oneindig en te "buitensporige" grenzen.
** Dit alles gebeurt zodanig dat alles keurig binnen de foutengrens
** blijft. CUP zal de kans aangeven van de afgekapte staart. BOUND krijgt
** voorlopig de waarde EPS, en zal zo aangepast worden dat uiteindelijk
** een goede schatting van de echt gemaakte fout gegeven wordt.
*/

if (eps<epsmin);
   errorlog "de opgegeven eps is te klein ";
   retp(2,2);
   end;
elseif (n<0 or n>7);
   errorlog "n heeft te extreme waarde";
   retp(2,2);
   end;
elseif not (((a .gt b) .or (inf .ne 2)) =3D=3D ones(n,1));
   errorlog "er zijn intervallen met slechte grenzen";
   retp(2,2);
   end;
endif;

bound=3Deps;
ept=3D(eps*0.15)/n;
z=3D-cdfninv(ept)+epsmin;
cup=3D1-cdfn(z);



/*
** BLADZIJDE C: inversie SIGMA en bepaling LRHPM
**
** Dit stuk inverteert SIGMA en de N-2 lower right hand principal
** minors. Alle inversies maken gebruik van INVPD. De LRHPM worden in
** een array van matrices gestopt (met naam LRHPM) via de VARPUT.
** Oproepen kunnen verder plaatsgrijpen via de VARGET. De controle
** variabele CNT vertelt of de VARPUT operatie geslaagd is. Buiten
** controle doeleinden heeft deze variabele geen enkel nut. Maar hem niet
** invoeren geeft telkens een output van 1 of 0 die geen betekenis heeft.
** Als N <=3D2 wordt dit stuk niet uitgevoerd. Aan het einde worden de
** variabelen gezet die betrekking hebben op het resterende bivariate
** geval.
*/


if (n>2);
   sinv=3Dinvpd(sigma);
   dett=3Ddetl;
   let lrhpm=3D lrhpm1 lrhpm2 lrhpm3 lrhpm4 lrhpm5 lrhpm6;
   tel=3D0;
   do while (tel<=3Dn-3);
      cnt=3Dvarput(invpd(sigma[(n-tel):n,(n-tel):n]),lrhpm[tel+1]);
      tel=3Dtel+1;
   endo;
   simp=3D1;
   if ((dett<0.05) or (eps<=3Depssim));simp=3D0;endif;
   prob=3D0;
   dett=3Dsinv[1,1]*sinv[2,2]-sinv[1,2]*sinv[1,2];
   sd[1]=3Dsqrt(sinv[2,2]/dett);
   sd[2]=3Dsqrt(sinv[1,1]/dett);
   rho=3D-sinv[2,1]/(sd[1]*sd[2]*dett);
   if (abs(rho)>1);
      errorlog "rho te groot";
      retp(2,2);
      end;
   endif;
endif;
nm2=3Dn-2;
nm1=3Dn-1;


/*
** BLADZIJDE D1: behandelen van de grenzen
**
** De grenzen worden afgekapt op [-Z,Z]. (Voor de betekenis van Z zie
** bladzijde B.) Dus A-waarden kleiner dan -Z worden vervangen door Z en
** analoog voor B-waarden. Waar een grens oneindig is wordt hetzelfde
** proc=82d=82 toegepast. We passen niet louter simpelweg de grenzen aan in
** de vectoren A en B maar schrijven de gecorrigeerde waarden naar C en D.
** Zo weten we waar wijzigingen zijn doorgevoerd en kunnen we het aantal
** gevallen waar dat gebeurd is tellen (T). De maximale fout per wijziging
** is EPLOS. Door grenzen aan te passen kan echter opnieuw D > C gaan
** optreden. Als dat gebeurd is wordt er gestopt. Tenslotte, indien N=3D1
** wordt dit op de triviale manier (met CDFN) opgelost. In dat geval
** wordt er dan ook gestopt.
*/


c=3Dmaxc(b'|((-z)*ones(1,n)));
d=3Dminc(a'|(( z)*ones(1,n)));

tel=3D1;
do while (tel<=3Dn);
   if (inf[tel]=3D=3D0); d[tel]=3Dz;  endif;
   if (inf[tel]=3D=3D1); c[tel]=3D-z; endif;
   tel=3Dtel+1;
endo;

t=3Dsumc((a .gt z) .or (inf .eq 0));
t=3Dt+sumc((b .lt -z) .or (inf .eq 1));
eplos=3Dt*cup;


if not (c lt d);
   errorlog "grenzen c en d incorrect ";
   retp(2,2);
   end;
endif;

if (n=3D=3D1);
   prob=3Dcdfn(d[1])-cdfn(c[1]);
   retp(prob,1E-12);
   end;
endif;


/*
** BLADZIJDE D2: behandelen van de grenzen (vervolg)
**
** We gaan de grenzen van de eerste en de tweede variabele aanpassen om er
** een resulterende bivariate verdeling van te maken. Zo N=3D2 wordt het
** eindresultaat hier onmiddellijk berekend. Als beide variabelen op ofwel
** -infty ofwel op +infty staan dan is het op te lossen met de procedure
** CDFBVN (of het complement ervan, dat wordt gegarandeerd door FAC)
** van GAUSS. Zo dat niet het geval is kan rectan de oplossing
** bieden. Door de bivariate verdelingsfunctie te gebruiken stijgt de
** nauwkeurigheid terug, daarom moet EPLOS gecorrigeerd worden. EPSI bevat
** dan de uiteindelijk toegestane fout die per variabele mag gemaakt
** worden.
*/


fac=3D1;
ipr=3D0;

if (inf[1]=3D=3D1 and inf[2]=3D=3D1);
   ipr=3D1;
   eplos=3Deplos-2*cup;
elseif (inf[1]=3D=3D0 and inf[2]=3D=3D0);
   fac=3D-1;
   ipr=3D1;
   d[1]=3Dc[1];
   d[2]=3Dc[2];
   eplos=3Deplos-2*cup;
endif;

if (n=3D=3D2);
   rho=3Dsigma[2,1];
   if (abs(rho)>1);
      errorlog "rho-bivariaat is te groot";
      retp(2,2);
      end;
   endif;

   if (ipr);
      y1=3Dd[1]*fac;
      y2=3Dd[2]*fac;
      prob=3Dcdfbvn(y1,y2,rho);
      retp(prob,bound);
      end;
   else;
      prob=3Drectan(c[1],d[1],c[2],d[2],rho);
      retp(prob,bound);
      end;
   endif;

endif;

epsi=3D(eps-eplos)/nm2;


/*
** BLADZIJDE E: regressiecoefficienten
**
** Regressiecoefficienten zoeken (BETA,BH). Grenzen zoeken voor de
** conditionele integralen (BINC). COND[L]=3Dconditionele variantie
** (CONDL[L] is de logaritme) van de variabele X(N-L+1), gegeven variabelen
** met grotere indices.
** HM bevat voor een bepaalde doorgang van L de hoofdminor die in
** LRHPM[L] gestockeerd wordt, dit om een gebruik van het toch wel
** plaatsrovende (niet tijdrovende) VARGET te vermijden.
*/


l=3D1;
do while (l<=3Dnm2);
   hm=3Dvarget(lrhpm[l]);
   cond[l,1]=3D1/sqrt(hm[1,1]);
   condl[l,1]=3Dln(cond[l]);

   i=3D1;
   do while (i<=3Dl);
      j=3D1;
      do while (j<=3Dl);
         beta[l,i]=3Dbeta[l,i]+sigma[n-l+j,n-l]*hm[j,l-i+1];
         j=3Dj+1;
      endo;
      i=3Di+1;
   endo;

   bh=3Dzeros(n-l,1);
   k=3Dn-l-1;
   bh[k+1]=3Dbeta[l,l];

   i=3D1;
   do while (i<=3Dk);
      j=3D1;
      do while (j<=3Dl);
         bh[i]=3Dbh[i]+sigma[j+n-l,n-l]*hm[j,1];
         j=3Dj+1;
      endo;
      i=3Di+1;
   endo;

   sigc=3D0;

   i=3D1;
   do while (i<=3Dl);
      j=3D1;
      do while (j<=3Di);
         if (i<=3Drows(bh));
            sigc=3Dsigc+bh[i]*bh[j]*sinv[i,j];
         endif;
         j=3Dj+1;
      endo;
      i=3Di+1;
   endo;

   binc[1,l]=3D1;
   binc[2,l]=3Dsqrt(sigc);
   binc[3,l]=3D2*sigc;
   binc[4,l]=3Dcons2*sigc*binc[2,l];
   binc[5,l]=3D12*sigc*sigc;
   if not (simp);
      binc[6,l]=3Dsigc*sigc*binc[1,l]*cons3;
      binc[7,l]=3Dsigc*sigc*sigc*cons4;
   endif;

   if (l=3D=3Dnm2);
      bh=3Dzeros(nm2,1);
      i=3D1;
      do while (i<=3Dnm2);
         j=3D1;
         do while (j<=3Dnm2);
            bh[i]=3Dbh[i]+sigma[j+2,1]*hm[l-i+1,j];
            j=3Dj+1;
         endo;
         i=3Di+1;
      endo;

   endif;

   l=3Dl+1;
endo;

l=3D1;


/*
** BLADZIJDE F: CO zal de wortels van de eerste 5, respectievelijk de
**              eerste 7 Hermite polynomen bevatten. Deze getallen zitten
**              al als constanten opgeslagen in de waarden CO en DDO.
*/

if simp;
   co=3Deo;
   iend=3D13;
   ien=3D5;
else;
   co=3Dddo;
   iend=3D25;
   ien=3D7;
endif;


/*
** BLADZIJDE G: initi=89le waarden.
**
** Nog enkele paramaters worden ge=8Bnitialiseerd. XSS bevat partieelsommen=
,
** gebruikt om de conditionele gemiddelden te berekenen.
*/


xss[.,1]=3Dzeros(7,1);
xm[1]=3D0;
prod[1]=3D1;
pr2[1]=3D1;

i=3D1;
do while (i<=3Dnm2);
   ni=3Dn-i+1;
   pr2[i+1]=3Dpr2[i]*(d[ni]-c[ni]);
   i=3Di+1;
endo;


/*
** BLADZIJDE H
**
** BINT[L] is een grens voor de fout op niveau L en dieper
** ANS[L]  is de geaccumuleerde integraal op niveau L
** PREP[L] bevat de integrand op niveau L
** PREB[L] begrenst de fout, geaccumuleerd op niveau L en dieper
*/


bint[nm1]=3D0;
lab90:
intvl[l]=3D2;
ans[l]=3D0;
bou4[l]=3D0;
bint[l]=3D0;
prep[l]=3D0;
preb[l]=3D0;

k=3D1;


/*
** BLADZIJDE I: Welke CO in het huidig interval ? S[L,.] bevat eindpunten
**              van de intervals waarop integrand op niveau L en afgeleiden
**              monotoon is. NUM[L] bevat het aantal intervallen.
**
** S[L,1] krijgt het eerste nulpunt. S[L,IEND+2] het laatste. NJS wordt
** gezet op de eerste CO[I]*COND[L] die binnen het aldus afgebakende interv=
al
** zit. MJS op de laatste. Als dat klaar is worden de tussenliggende punten
** ook opgevuld zodat gans S[L,.] opgevuld is van 1 tot NUML. Op het ogenbl=
ik
** dat het aantal intervals 2 is (NUM[L]=3D=3D2) wordt er uiteraard gestopt=
.
*/


nl=3Dn-l+1;
s[l,1]=3Dc[nl]-xm[l];
s[l,iend+2]=3Dd[nl]-xm[l];
num[l]=3Diend+2;

i=3D1;
ok=3D0;
do while ((i<=3Diend) and not(ok));
   njs=3Di;
   if (s[l,1]<co[i]*cond[l]);
      ok=3D1;
   else;
      num[l]=3Dnum[l]-1;
   endif;
   i=3Di+1;
endo;

if (num[l]/=3D2);

   i=3Dnjs;
   ok=3D0;
   do while ((i<=3Diend) and (not ok));
      mjs=3Diend-i+njs;
      if (s[l,iend+2]>=3Dco[mjs]*cond[l]);
         ok=3D1;
      else;
         num[l]=3Dnum[l]-1;
      endif;
      i=3Di+1;
   endo;

   if (num[l]/=3D2);
      i=3Dnjs;
      do while (i<=3Dmjs);
         inj=3Di-njs+2;
         s[l,inj]=3Dco[i]*cond[l];
         i=3Di+1;
      endo;
   endif;

endif;

numl=3Dnum[l];
s[l,numl]=3Ds[l,iend+2];


/*
** BLADZIJDE J
**
** EP[L]    bovengrens voor de toelaatbare fout op niveau L
** R[K,L]   rechts eindpunt huidig deelinterval
** FE[L]    links eindpunt huidig deelinterval
** IND[L]-1 duidt aan welk punt van Newton-Cotes we behandelen
**
** Begrenzen van de afgeleiden aan links eindpunt huidig interval :
** BL[1,L]  grens i-de afgeleide normale dichtheid op niveau L aan links
**          eindpunt.
** Analoge bewerkingen voor het rechtse eindpunt gebeuren op de volgende
** bladzijde.
*/


ep[l]=3Depsi/pr2[l+1];
r[1,l]=3Ds[l,2];
ind[l]=3D6;

fe[l]=3Ds[l,1];
t=3Dfe[l]/cond[l];
bl[1,l]=3Dphi(t,condl[l]);
bl[2,l]=3Dbl[1,l]*abs(t/cond[l]);
bl[3,l]=3Dbl[1,l]*f2(t,cond[l]);
bl[4,l]=3Dbl[1,l]*f3(t,cond[l]);
bl[5,l]=3Dbl[1,l]*f4(t,cond[l]);

if not (simp);
   bl[6,l]=3Dbl[1,l]*f5(t,cond[l]);
   bl[7,l]=3Dbl[1,l]*f6(t,cond[l]);
endif;


/*
** BLADZIJDE K: Begrenzen van de afgeleiden aan het rechts eindpunt van het
**              huidig interval.
**
** BR[I][K,L] is analoog aan BL[I,L], nu is er echter wel een K-index.
** Er moet hier een technisch probleempje overwonnen worden. GAUSS laat nie=
t
** toe van met meerdimensionale matrices te werken. Daarom moeten we nu een
** array construeren van namen van matrices.
** Dat is gebeurd in de initialisatie.
*/


lab100:
t=3Dr[k,l]/cond[l];
getal=3Dphi(t,condl[l]);

hulp=3Dvarget(br[1]);
hulp[k,l]=3Dgetal;
cnt=3Dvarput(hulp,br[1]);

hulp=3Dvarget(br[2]);
hulp[k,l]=3Dgetal*abs(t/cond[l]);
cnt=3Dvarput(hulp,br[2]);

hulp=3Dvarget(br[3]);
hulp[k,l]=3Dgetal*f2(t,cond[l]);
cnt=3Dvarput(hulp,br[3]);

hulp=3Dvarget(br[4]);
hulp[k,l]=3Dgetal*f3(t,cond[l]);
cnt=3Dvarput(hulp,br[4]);

hulp=3Dvarget(br[5]);
hulp[k,l]=3Dgetal*f4(t,cond[l]);
cnt=3Dvarput(hulp,br[5]);

if not (simp);

   hulp=3Dvarget(br[6]);
   hulp[k,l]=3Dgetal*f5(t,cond[l]);
   cnt=3Dvarput(hulp,br[6]);

   hulp=3Dvarget(br[7]);
   hulp[k,l]=3Dgetal*f6(t,cond[l]);
   cnt=3Dvarput(hulp,br[7]);

endif;

lab104:
r[k+1,l]=3D(fe[l]+r[k,l])*0.5;
bou5[l]=3Dep[l]*(r[k,l]-s[l,1]);
del[2,l]=3Dr[k+1,l]-fe[l];


/*
** BLADZIJDE L:  testen van de grenzen voor de Newton-Cotes regels
**
** We zijn bezig te testen voor de l-de variabele welke integratieformule
** gebruikt wordt. De drie regels worden getest (progressief). Als er =82=
=82n
** goed is stopt de test (er wordt naar LAB200 gewipt). Voor Simpson is
** er een alternatief: als de variabele SIMP true is (waarde 1) komt de
** Boole-regel niet in aanmerking (dit heeft te maken met de determinant
** van SIGMA). In dat geval voorziet men het overslaan van Booles regel
** en wordt er naar LAB130 gewipt. Als alle testen falen zal het interval
** verder opgesplitst worden en springen we naar LAB100.
** Vanaf nu wordt er vrij intensief gebruik gemaakt van GOTO statements.
** Er werd naar gestreefd enerzijds de leesbaarheid en de overzichtelijkhei=
d
** maximaal te houden en anderzijds de implementatie van Shervish in Fortra=
n
** zo dicht mogelijk te blijven benaderen. We hebben geoordeeld dat het
** programma hoe dan ook als gestructureerd dient beschouwd.
*/


/* test voor de trapeziumregel */

del[3,l]=3D2*del[2,l];
hulp1=3Dvarget(br[1]);
hulp2=3Dvarget(br[2]);
hulp3=3Dvarget(br[3]);
bou1=3Dmaxc(hulp1[k,l]|bl[1,l])*binc[3,l];
bou1=3Dbou1+2*maxc(hulp2[k,l]|bl[2,l])*binc[2,l];
bou1=3Dbou1+maxc(hulp3[k,l]|bl[3,l])*binc[1,l];
bou3=3Dbou4[l]+bou1*(del[3,l]^3)*prod[l]/12;
itype[l]=3D3;

if (bou3<=3Dbou5[l]); goto lab200; endif;

/* test voor de Simpsonregel, indien nodig */

bou1=3D0;
ij=3D1;
do while (ij<=3D5);
   jk=3D6-ij;
   hulp=3Dvarget(br[ij]);
   bou2=3Dmaxc(hulp[k,l]|bl[ij,l]);
   bou1=3Dbou1+bou2*binc[jk,l]*bcs[ij];
   ij=3Dij+1;
endo;
bou3=3Dbou4[l]+bou1*(del[2,l]^5)*prod[l]/90;
itype[l]=3D2;

if (bou3<=3Dbou5[l]); goto lab200; endif;

if simp; goto lab130; endif;

/* test voor de Booleregel, indien nodig */

del[1,l]=3D0.5*del[2,l];
bou1=3D0;
itype[l]=3D1;
ij=3D1;
do while (ij<=3D7);
   jk=3D8-ij;
   hulp=3Dvarget(br[ij]);
   bou2=3Dmaxc(hulp[k,l]|bl[ij,l]);
   bou1=3Dbou1+bou2*binc[jk,l]*bcn[ij];
   ij=3Dij+1;
endo;
bou3=3Dbou4[l]+bou1*(del[1,l]^7)*prod[l]*8/945;

if (bou3<=3Dbou5[l]); goto lab200; endif;


/*
** BLADZIJDE M: het interval verder onderverdelen als alle testen faalden
**
** Hier komen we aan als de drie regels een negatief antwoord geven. We
** gaan dan het interval verder onderverdelen en opnieuw de testen uitvoere=
n.
** Hiertoe dient K vergroot. Uit veiligheidsoverwegingen is er een stop
** voorzien als K te groot. In dat geval wordt ook een foutmelding gemaakt.
** Indien inderdaad de K dient vergroot, dan wordt er naar LAB100 gespronge=
n.
** Merk op dat men op deze plaats niet passeert als =82=82n van de testen
** "geslaagd" is. In dat geval zijn we al naar LAB200 gewipt.
*/


lab130:
k=3Dk+1;
if (k>35);
   errorlog "de k-waarde is te groot geworden";
   retp(2,2);
   end;
endif;
goto lab100;


/*
** BLADZIJDE N: vertrekken naar juiste plaats volgens IND en ITYPE
*/

lab200:
bint[l]=3Dbint[l]+bou3-bou4[l];
bou4[l]=3Dbou3;
ksa[l]=3Dk;
if (ind[l]=3D=3D6); goto lab202; endif;
if (itype[l]=3D=3D1);
   goto lab205;
elseif (itype[l]=3D=3D2);
   goto lab206;
elseif (itype[l]=3D=3D3);
   goto lab210;
endif;


/*
** BLADZIJDE P: bepaling XS, naar gelang van IND[] en ITYPE.
**
** Indien mogelijk wordt een niveau dieper gewipt om daar de integratie
** verder te zetten. Zoniet zal de uiteindelijke evaluatie begonnen worden.
*/


lab202:
ind[l]=3D5;
xs=3Dfe[l];
fact[l]=3Dbl[1,l];
lab203:
xss[nm1,l+1]=3Dxss[nm1,l]+bh[l]*(xs+xm[l]);
ll=3Dl;
do while (ll<=3Dnm2);
   xss[ll,l+1]=3Dxss[ll,l]+beta[ll,l]*(xs+xm[l]);
   ll=3Dll+1;
endo;
if (l=3D=3Dnm2); goto lab300; endif;
xm[l+1]=3Dxss[l,l+1];
prod[l+1]=3Dprod[l]*fact[l];
l=3Dl+1;
goto lab90;

lab205:
ind[l]=3D4;
k=3Dksa[l];
xs=3D0.5*(fe[l]+r[k+1,l]);
goto lab207;

lab206:
ind[l]=3D3;
k=3Dksa[l];
xs=3Dr[k+1,l];

lab207:
t=3Dxs/cond[l];
fact[l]=3Dphi(t,condl[l]);
goto lab203;

lab208:
ind[l]=3D2;
k=3Dksa[l];
xs=3D0.5*(r[k,l]+r[k+1,l]);
goto lab207;

lab210:
ind[l]=3D1;
k=3Dksa[l];
xs=3Dr[k,l];
hulp=3Dvarget(br[1]);
fact[l]=3Dhulp[k,l];
goto lab203;


/*
**    BLADZIJDE Q: Evaluatie conditionele bivariate integraal op diepste
**                 niveau. Als L=3D1 wordt hier gestopt.
*/


lab300:
x1=3D-fac*(xss[nm1,nm1]-d[1])/sd[1];
x2=3D-fac*(xss[nm2,nm1]-d[2])/sd[2];
l=3Dnm1;
y1=3D-(xss[nm1,nm1]-c[1])/sd[1];
y2=3D-(xss[nm2,nm1]-c[2])/sd[2];

if (ipr);
   ans[l]=3Dcdfbvn(x1,x2,rho);
else;
   ans[l]=3Drectan(x1,y1,x2,y2,rho);
endif;


lab310:
if (l=3D=3D1);
   prob=3Dans[1];
   bound=3Dbint[1]+eplos;
   retp(prob,bound);
   end;
endif;
l=3Dl-1;


/*
** BLADZIJDE R: integratie op huidig niveau verderzetten
*/


indl=3Dind[l];
numl=3Dnum[l];
ity=3Ditype[l];
temp=3Dfact[l]*ans[l+1];
temb=3Dbint[l+1];
fsa=3D1;

if (indl=3D=3D1);
   tem=3Dtemp;
   temp=3Dtemp+prep[l];
   temb=3Dtemb+preb[l];
   prep[l]=3Dtem;
   preb[l]=3Dbint[l+1];
   fsa=3D2;
endif;

ans[l]=3Dans[l]+coeff[indl,ity]*temp*del[ity,l];
bint[l]=3Dbint[l]+coeff[indl,ity]*temb*del[ity,l];


/*
** BLADZIJDE S
**
** Gebruik maken van de niet geaccumuleerde fout op niveau L+1.
** De onderverdeling ongedaan maken op niveau L.
*/


ep[l]=3Dep[l]+del[ity,l]*coeff[indl,ity]*(fsa*(nm2-l)*epsi/pr2[l+1]-temb)/
      (s[l,numl]-s[l,1]);

if (indl/=3D1);
   igo=3Dindl-(1+(itype[l]*(itype[l]-1))/2);
   if igo=3D=3D1;
      goto lab210;
   elseif igo=3D=3D2;
      goto lab208;
   elseif igo=3D=3D3;
      goto lab206;
   elseif igo=3D=3D4;
      goto lab205;
   endif;
endif;

k=3Dksa[l];

i=3D1;
do while (i<=3Dien);
   hulp=3Dvarget(br[i]);
   bl[i,l]=3Dhulp[k,l];
   i=3Di+1;
endo;

ind[l]=3D5;
fe[l]=3Dr[k,l];


if (k/=3D1);
   k=3Dk-1;
   goto lab104;
endif;

if (intvl[l]=3D=3Dnum[l]);
   goto lab310;
endif;

intvl[l]=3Dintvl[l]+1;
intl=3Dintvl[l];
r[1,l]=3Ds[1,intl];

goto lab100;

endp;


