[6557] in s-news-athena
Re: Cumulative multivariate normal distribution.
daemon@ATHENA.MIT.EDU (Huiman Xie Barnhart)
Fri Mar 3 10:00:25 1995
Date: Fri, 3 Mar 1995 08:51:37 -0500 (EST)
From: Huiman Xie Barnhart <hxb@panda.sph.emory.edu>
To: s-news@utstat.toronto.edu
Cc: "Byers, Bob" <RHB1@CIDHIV1.EM.CDC.GOV>
Lots of thanks to all who send me information and programs which
calculate the cumulative multivariate normal distribution.
See below for all the responses I got.
Huiman
--------------------------------------------
From mtan@bio.ri.ccf.org
Date: Wed, 1 Mar 1995 17:18:15 -0500
From: Ming Tan <mtan@bio.ri.ccf.org>
To: hxb@panda.sph.emory.edu
Thanks to many people who sent me the bivariate cdf functions.
Here are the three functions to calculate bivariate prob :
(i) Splus calling C
(ii) Splus calling Fortran
(iii) pure Splus
Itried all of they three. All of them worked.
Some system related changes have to be made for the C-code.
Ming Tan, Ph.D.
Department of Biostatistics and Epidemiology
The Cleveland Clinic Foundation
----------------------
The FOLLOWING is the Splus calling C
>From tim@fnmcps.fandm.edu Sun Mar 7 09:13:32 1993
Date: 07 Mar 1993 09:09:43 -0500
From: "Prof. Tim Hesterberg" <tim@fnmcps.fandm.edu>
Subject: bivariate normal prob
To: mtan@bio.ri.ccf.org
Content-Transfer-Encoding: 7BIT
Content-Length: 7733
X-Lines: 334
Status: RO
This message contains two files, which are the S and C code for
bivariate normal probaibilities. You will have to figure out the
correct compilation for ANSI C for your machine.
I use S-Plus.
The C code calls a function "not.loaded()", which is a modified version of
the S-Plus function "is.loaded()". Here is "not.loaded()":
not.loaded <- function(symbol)
.C("S_get_entry",
as.character(symbol.C(symbol)),
integer(1))[[2]] <= 0
First, my file "binorm2.S", which contains code for function
pbinorm2()
"pbinorm2"<- function(x, y, rho, mu.x=0, mu.y=0, sigma.x=1, sigma.y=1 )
{
## Returns bivariate normal probability. Input may be vectors.
if( not.loaded( "pbinorm2" ))
dyn.load("/users/tim/Scode/binorm2.o")
if( any( abs(rho) > 1 )) stop("correlations must be between -1 and 1")
if( !missing(mu.x) || !missing(mu.y)) x _ (x-mu.x)/sigma.x
if( !missing(mu.y) || !missing(mu.y)) y _ (y-mu.y)/sigma.y
l _ max( length(x), length(y), length(rho) )
.C( "pbinorm2", as.integer(l),
x = as.double( rep( x,length=l) ),
y = as.double( rep( y,length=l) ),
rho = as.double( rep( rho,length=l) ),
results = double( l ))$results
}
Now, my file "binorm2.c", which includes the code from statlib, and a
small interface routine I wrote at the bottom.
/* From statlib, file general/gaut.c, by Ajay Shah */
/* Appears to be ANSI C; compile using occ -A -c binorm2.c */
/* Interface routine "pbinorm2" to S by Tim Hesterberg is at bottom. */
#include "math.h"
/* Frequently used numerical constants: */
#define OneUponSqrt2Pi .39894228040143267794
#define twopi 6.283195307179587
#define LnSqrt2Pi -0.9189385332046727417803296
#define SQRT2 1.414213562373095049
#define SQRTPI 1.772453850905516027
/* ---------------------------------------------------------------------------
UNIVARIATE NORMAL PROBABILITY
---------------------------------------------------------------------------*/
#define UPPERLIMIT 20.0
/* I won't return either of univariate normal density or
probability when x < -UPPERLIMIT or x > UPPERLIMIT. */
#define P10 242.66795523053175
#define P11 21.979261618294152
#define P12 6.9963834886191355
#define P13 -.035609843701815385
#define Q10 215.05887586986120
#define Q11 91.164905404514901
#define Q12 15.082797630407787
#define Q13 1.0
#define P20 300.4592610201616005
#define P21 451.9189537118729422
#define P22 339.3208167343436870
#define P23 152.9892850469404039
#define P24 43.16222722205673530
#define P25 7.211758250883093659
#define P26 .5641955174789739711
#define P27 -.0000001368648573827167067
#define Q20 300.4592609569832933
#define Q21 790.9509253278980272
#define Q22 931.3540948506096211
#define Q23 638.9802644656311665
#define Q24 277.5854447439876434
#define Q25 77.00015293522947295
#define Q26 12.78272731962942351
#define Q27 1.0
#define P30 -.00299610707703542174
#define P31 -.0494730910623250734
#define P32 -.226956593539686930
#define P33 -.278661308609647788
#define P34 -.0223192459734184686
#define Q30 .0106209230528467918
#define Q31 .191308926107829841
#define Q32 1.05167510706793207
#define Q33 1.98733201817135256
#define Q34 1.0
double pnorm1(double x)
{
int sn;
double R1, R2, R3, y, y2, y3, y4, y5, y6, y7;
double erf, erfc, z, z2, z3, z4;
double phi;
if (x < -UPPERLIMIT) return 0.0;
if (x > UPPERLIMIT) return 1.0;
y = x / SQRT2;
if (y < 0) {
y = -y;
sn = -1;
}
else
sn = 1;
y2 = y * y;
y4 = y2 * y2;
y6 = y4 * y2;
if(y < 0.46875) {
R1 = P10 + P11 * y2 + P12 * y4 + P13 * y6;
R2 = Q10 + Q11 * y2 + Q12 * y4 + Q13 * y6;
erf = y * R1 / R2;
if (sn == 1)
phi = 0.5 + 0.5*erf;
else
phi = 0.5 - 0.5*erf;
}
else
if (y < 4.0) {
y3 = y2 * y;
y5 = y4 * y;
y7 = y6 * y;
R1 = P20 + P21 * y + P22 * y2 + P23 * y3 +
P24 * y4 + P25 * y5 + P26 * y6 + P27 * y7;
R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 +
Q24 * y4 + Q25 * y5 + Q26 * y6 + Q27 * y7;
erfc = exp(-y2) * R1 / R2;
if (sn == 1)
phi = 1.0 - 0.5*erfc;
else
phi = 0.5*erfc;
}
else {
z = y4;
z2 = z * z;
z3 = z2 * z;
z4 = z2 * z2;
R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4;
R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4;
erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2));
if (sn == 1)
phi = 1.0 - 0.5*erfc;
else
phi = 0.5*erfc;
}
return phi;
}
/* ---------------------------------------------------------------------------
UNIVARIATE NORMAL DENSITY
---------------------------------------------------------------------------*/
double dnorm1(double x)
{
if (x < -UPPERLIMIT) return 0.0;
if (x > UPPERLIMIT) return 0.0;
return OneUponSqrt2Pi*exp(-0.5*x*x);
}
/* ---------------------------------------------------------------------------
LN OF UNIVARIATE NORMAL DENSITY
---------------------------------------------------------------------------*/
double lndnorm1(double x)
{
return LnSqrt2Pi - (0.5*x*x);
}
/*---------------------------------------------------------------------------
BIVARIATE NORMAL PROBABILITY
---------------------------------------------------------------------------*/
#define con (twopi / 2.0) * 10.0e-10
double bivnor(double ah, double ak, double r)
{
/*
based on alg 4628 comm. acm oct 73
gives the probability that a bivariate normal exceeds (ah,ak).
gh and gk are .5 times the right tail areas of ah, ak under a n(0,1)
Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand.
*/
double a2, ap, b, cn, conex, ex, g2, gh, gk, gw, h2, h4, rr, s1, s2,
sgn, sn, sp, sqr, t, temp, w2, wh, wk;
int is;
temp = -ah;
gh = pnorm1(temp);
gh = gh / 2.0;
temp = -ak;
gk = pnorm1(temp);
gk = gk / 2.0;
b = 0;
if (r==0)
b = 4*gh*gk;
else {
rr = 1-r*r;
if (rr<0)
return;
if (rr!=0) {
sqr = sqrt(rr);
if (ah!=0) {
b = gh;
if (ah*ak<0)
b = b-.5;
else if (ah*ak==0)
goto label10;
}
else if (ak==0) {
b = atan(r/sqr)/twopi+.25;
goto label50;
}
b = b+gk;
if (ah==0)
goto label20;
label10:
wh = -ah;
wk = (ak/ah-r)/sqr;
gw = 2*gh;
is = -1;
goto label30;
label20:
do {
wh = -ak;
wk = (ah/ak-r)/sqr;
gw = 2*gk;
is = 1;
label30:
sgn = -1;
t = 0;
if (wk!=0) {
if (fabs(wk)>=1)
if (fabs(wk)==1) {
t = wk*gw*(1-gw)/2;
goto label40;
}
else {
sgn = -sgn;
wh = wh*wk;
g2 = pnorm1(wh);
wk = 1/wk;
if (wk<0)
b = b+.5;
b = b-(gw+g2)/2+gw*g2;
}
h2 = wh*wh;
a2 = wk*wk;
h4 = h2*.5;
ex = 0;
if (h4<150.0)
ex = exp(-h4);
w2 = h4*ex;
ap = 1;
s2 = ap-ex;
sp = ap;
s1 = 0;
sn = s1;
conex = fabs(con/wk);
do {
cn = ap*s2/(sn+sp);
s1 = s1+cn;
if (fabs(cn)<=conex)
break;
sn = sp;
sp = sp+1;
s2 = s2-w2;
w2 = w2*h4/sp;
ap = -ap*a2;
} while (1);
t = (atan(wk)-wk*s1)/twopi;
label40:
b = b+sgn*t;
}
if (is>=0)
break;
} while(ak!=0);
}
else if (r>=0)
if (ah>=ak)
b = 2*gh;
else
b = 2*gk;
else if (ah+ak<0)
b = 2*(gh+gk)-1;
}
label50:
if (b<0)
b = 0;
if (b>1)
b = 1;
return(b);
}
void pbinorm2( long *m, double *x, double *y, double *r, double *result )
{
/* Call bivnor *m times, using corresponding values in x, y, r. */
/* Return results in *result. */
/* This gives the cumulative distribution function, hence the "-"'s. */
int i;
for(i=0;i< *m;i++) *result++ = bivnor( -*x++, -*y++, *r++);
}
Hope this helps.
Tim Hesterberg, tim@fnmcps.fandm.edu, (717)291-3923
Mathematics Department, Franklin & Marshall College, Lancaster, PA 17604-3003
THE FOLLOWING is the function calling FORTRAN from Daniel Heitjan:
>From dheitjan@biostats.hmc.psu.edu Thu Mar 4 17:06:35 1993
Date: Thu, 4 Mar 93 17:07:54 EST
From: dheitjan@biostats.hmc.psu.edu (Daniel F. Heitjan)
To: mtan@bio.ri.ccf.org
Subject: Re: bivariate normal integral
Content-Length: 500
Status: RO
X-Lines: 17
Here is the S. I'll send the fortran under
separate cover.
fun.bvnprob_function(mu,sig,ah,ak) {
## Call the Fortran program that computes the BVN area. This
## version computes the area in the northeast rectangle.
## Daniel F. Heitjan, 23 February 1991
ah_as.double(ah)
ak_as.double(ak)
mu_c(mu)
xmu1_as.double(mu[1])
xmu2_as.double(mu[2])
v1_as.double(sig[1,1])
v2_as.double(sig[2,2])
c12_as.double(sig[1,2])
prob_.Fortran('bvnprb',ah,ak,xmu1,xmu2,v1,v2,c12,as.double(0))
prob[[8]] }
>From dheitjan@biostats.hmc.psu.edu Thu Mar 4 17:09:40 1993
Date: Thu, 4 Mar 93 17:11:09 EST
From: dheitjan@biostats.hmc.psu.edu (Daniel F. Heitjan)
To: mtan@bio.ri.ccf.org
Subject: Re: bivariate normal integral--fortran parts
Content-Length: 5879
Status: RO
X-Lines: 228
subroutine bvnprb(ah,ak,xmu1,xmu2,v1,v2,c12,prob)
c compute the probability that a bivariate normal with mean vectorc
(xmu1,xmu2), variances v1 and v2 and covariance c12 has x1>ah
c and x2>ak.
c daniel f. heitjan, 17 march 1991
implicit real*8 (a-h,o-z)
ahi=(ah-xmu1)/v1**0.5
aki=(ak-xmu2)/v2**0.5
r=c12/(v1*v2)**0.5
prob=bivnor(ahi,aki,r)
return
end
subroutine calerf(arg,result,jint)
integer i,jint
real*8 a,arg,b,c,d,four,half,p,one,q,result,sqrpi,
# two,thresh,x,xmax,xden,xnum,xsmall,y,ysq,zero
dimension a(5),b(4),c(9),d(8),p(6),q(5)
c mathematical constants
data four,one,half,two,zero/4.0d0,1.0d0,0.5d0,2.0d0,0.0d0/
data sqrpi/5.6418958354775628695d-1/,thresh/0.46875d0/
c machine-dependent parameters
data xsmall/4.2d-16/,xmax/9.269d0/
c coefficients for approximation to derf in first interval
data a/3.16112374387056560d00,1.13864154151050156d02,
# 3.77485237685302021d02,3.20937758913846947d03,
# 1.85777706184603153d-1/
data b/2.36012909523441209d01,2.44024637934444173d02,
# 1.28261652607737228d03,2.84423683343917062d03/
c coefficients for approximation to derfc in second interval
data c/5.64188496988670089d-1,8.88314979438837594d0,
# 6.61191906371416295d01,2.98635138197400131d02,
# 8.81952221241769090d02,1.71204761263407058d03,
# 2.05107837782607147d03,1.23033935479799725d03,
# 2.15311535474403846d-8/
data d/1.57449261107098347d01,1.17693950891312499d02,
# 5.37181101862009858d02,1.62138957456669019d03,
# 3.29079923573345963d03,4.36261909014324716d03,
# 3.43936767414372164d03,1.23033935480374942d03/
c coefficients for approximation to derfc in third interval
data p/3.05326634961232344d-1,3.60344899949804439d-1,
# 1.25781726111229246d-1,1.60837851487422766d-2,
# 6.58749161529837803d-4,1.63153871373020978d-2/
data q/2.56852019228982242d00,1.87295284992346047d00,
# 5.27905102951428412d-1,6.05183413124413191d-2,
# 2.33520497626869185d-3/
x=arg
y=dabs(x)
if (y.gt.four) go to 200
if (y.gt.thresh) go to 100
c evaluate erf for abs(x).le.0.46875
ysq=zero
if (y.gt.xsmall) ysq=y*y
xnum=a(5)*ysq
xden=ysq
do 20 i=1,3
xnum=(xnum+a(i))*ysq
xden=(xden+b(i))*ysq
20 continue
result=x*(xnum+a(4))/(xden+b(4))
if (jint.ne.0) result=one-result
go to 800
c evaluate erfc for 0.46875.lt.abs(x).le.4.0
100 ysq=y*y
xnum=c(9)*y
xden=y
do 120 i=1,7
xnum=(xnum+c(i))*y
xden=(xden+d(i))*y
120 continue
result=dexp(-ysq)*(xnum+c(8))/(xden+d(8))
go to 300
c evaluate erfc for abs(x).gt.4.0
200 result=zero
if (y.ge.xmax) go to 300
220 ysq=one/(y*y)
xnum=p(6)*ysq
xden=ysq
do 240 i=1,4
xnum=(xnum+p(i))*ysq
xden=(xden+q(i))*ysq
240 continue
result=ysq*(xnum+p(5))/(xden+q(5))
result=(dexp(-y*y)/y)*(sqrpi-result)
c fix up for neg. arg.,erf,etc.
300 if (jint.eq.0) go to 350
if (x.lt.zero) result=two-result
go to 800
350 result=(half-result)+half
if (x.lt.zero) result=-result
800 return
end
real*8 function derf(x)
c program to compute the error function
c author-w. j. cody
c date-january 8,1985
integer jint
real*8 x,result
jint=0
call calerf(x,result,jint)
derf=result
return
end
double precision function bivnor(ah,ak,r)
double precision ah,ak,r
double precision twopi,b,xah,xak,gh,gk,rr,gauss,derf,h2,
+ a2,h4,ex,w2,ap,s2,sp,s1,sn,sgn,sqr,con,wh,wk,gw,
+ t,g2,conex,cn,two,zero,one,four,quart,half,explim
integer idig,is
data two/2.d0/,zero/0.d0/,one/1.d0/,four/4.d0/,quart/0.25d0/,
+ half/0.5d0/
gauss(t)=(one+derf(t/sqrt(two)))/two
data twopi/6.2831 85307 17958 7d0/,idig/15/,explim/80.d0/
b=zero
if (ah.eq.zero) then
xah=ak
xak=ah
else
xah=ah
xak=ak
end if
gh=gauss(-xah)/two
gk=gauss(-xak)/two
if (r.eq.zero) then
b=four*gh*gk
go to 350
end if
rr=one-r*r
if (rr.lt.zero) then
c write(*, *)'error in bivnor, r=', r
go to 390
end if
if (rr.gt.zero) go to 100
if (r.ge.zero) go to 70
if (xah+xak.ge.zero) go to 350
b=two*(gh+gk)-one
go to 350
70 if (xah-xak.lt.zero) then
b=two*gk
else
b=two*gh
end if
go to 350
100 sqr=sqrt(rr)
if (idig.eq.15) then
con=twopi*1.d-15/two
else
con=twopi/two/10**idig
end if
if (xah.ne.zero) go to 170
if (xak.ne.zero) go to 190
b=atan(r/sqr)/twopi+quart
go to 350
170 b=gh
if (xah*xak) 180, 200, 190
180 b=b-half
190 b=b+gk
200 wh=-xah
wk=(xak/xah-r)/sqr
gw=two*gh
is=-1
210 sgn=-one
t=zero
if (wk.eq.zero) go to 320
if (abs(wk)-one) 270, 230, 240
230 t=wk*gw*(one-gw)/two
go to 310
240 sgn=-sgn
wh=wh*wk
g2=gauss(wh)
wk=one/wk
if (wk.lt.zero) b=b+half
b=b-(gw+g2)/two+gw*g2
270 h2=wh*wh
a2=wk*wk
h4=h2/two
if (h4.lt.explim) then
ex=exp(-h4)
else
ex=zero
end if
w2=h4*ex
ap=one
s2=ap-ex
sp=ap
s1=zero
sn=s1
conex=abs(con/wk)
go to 290
280 sn=sp
sp=sp+one
s2=s2-w2
w2=w2*h4/sp
ap=-ap*a2
290 cn=ap*s2/(sn+sp)
s1=s1+cn
if (abs(cn)-conex.gt.zero) go to 280
t=(atan(wk)-wk*s1)/twopi
310 b=b+sgn*t
320 if (is.ge.0) go to 350
if (xak.ne.zero) then
wh=-xak
wk=(xah/xak-r)/sqr
gw=two*gk
is=1
go to 210
end if
350 if (b.lt.zero) b=zero
if (b.gt.one) b=one
390 bivnor=b
return
end
THE FOLLOWING is the pure splus function sent by
Prof. John Kershaw <kershaw@math.unb.ca> sent S code.
P.bi.Norm<-function(Zx,Zy,rho,conv=.000001)
{
if(abs(rho) > 1){stop("| rho | > 1")}
Bprob_0
if(Zx==0)
{
tmp_Zx
Zx_Zy
Zy_tmp
}
UPx_pnorm(-Zx)/2
UPy_pnorm(-Zy)/2
if(rho==0)
{
Bprob_4*UPx*UPy
}
else
{
RR_1-(rho^2)
if(RR==0)
{
if(rho > 0)
{
if((Zx-Zy) < 0)
{
Bprob_2*UPy
}
else
{
Bprob_2*UPx
}
}
else
{
if((Zx+Zy) >= 0)
{
Bprob_0
}
else
{
Bprob_(2*(UPx+UPy))-1
}
}
}
else
{
CONT_F
conn_pi*conv
sqrr_sqrt(RR)
if(Zx != 0)
{
Bprob_UPx
if((Zx*Zy) != 0)
{
if((Zx*Zy) < 0){Bprob_Bprob-.5}
Bprob_Bprob+UPy
}
}
else
{
if(Zy != 0)
{
Bprob_Bprob+UPy
}
else
{
Bprob_(atan(rho/sqrr)/(2*pi))+.25
CONT_T
}
}
wh_-Zx
wk_((Zy/Zx)-rho)/sqrr
gw_2*UPx
swt_F
while(!CONT)
{
SGN_-1
TTT_0
if(wk==0)
{
if(swt){break}
if(Zy==0){break}
wh_-Zy
wk_((Zx/Zy)-rho)/sqrr
gw_2*UPy
swt_T
}
if(wk==0){break}
tst_abs(wk)-1
if(tst==0)
{
TTT_wk*gw*(1-gw)/2
Bprob_Bprob+(SGN*TTT)
if(swt){break}
if(Zy==0){break}
wh_-Zy
wk_((Zx/Zy)-rho)/sqrr
gw_2*UPy
swt_T
}
else
{
if(tst > 0)
{
SGN_-SGN
wh_wh*wk
g2_pnorm(wh)
wk_1/wk
if(wk < 0){Bprob_Bprob+.5}
Bprob_Bprob-((gw+g2)/2)+(gw*g2)
}
h2_wh*wh
a2_wk*wk
h4_h2/2
if(h4 < 80){EX_exp(-h4)}else{EX_0}
w2_h4*EX
AP_1
s2_AP-EX
SP_1
s1_0
SN_s1
conex_abs(conn/wk)
CN_AP*s2/(SN+SP)
s1_s1+CN
while((abs(CN)-conex) > conv)
{
SN_SP
SP_SP+1
s2_s2-w2
w2_w2*h4/SP
AP_-AP*a2
CN_AP*s2/(SN+SP)
s1_s1+CN
}
TTT_(atan(wk)-(wk*s1))/(2*pi)
Bprob_Bprob+(SGN*TTT)
if(swt){break}
if(Zy==0){break}
wh_-Zy
wk_((Zx/Zy)-rho)/sqrr
gw_2*UPy
swt_T
}
}
}
}
if(Bprob < 0){Bprob_0}
if(Bprob > 1){Bprob_1}
return(Bprob)
}
pbinorm.region<-function(x1,x2,y1,y2,rho)
{
P1_P.bi.Norm(x1,y1,rho)
P2_P.bi.Norm(x1,y2,rho)
P3_P.bi.Norm(x2,y1,rho)
P4_P.bi.Norm(x2,y2,rho)
PPP_P1-P2-P3+P4
return(PPP)
}
I timed the three functions :
> unix.time(bivnorms(1.2,1.8,0.2)) # pure Splus
[1] 0.13333368 0.08333302 0.00000000 0.00000000 0.00000000
> unix.time(bivnormf(1.2,1.8,0.2)) # Splus calling C
[1] 0.09999990 0.01666641 0.00000000 0.00000000 0.00000000
> unix.time(bivnormc(1.2,1.8,0.2)) #Splus calling Fortran
[1] 0.01666641 0.00000000 0.00000000 0.00000000 0.00000000
Ming Tan
The Cleveland Clinic
mtan@bio.ri.ccf.org
From acbsjack@gosnell.spc.uchicago.edu
Date: Wed, 1 Mar 1995 18:55:44 -0600 (CST)
From: Simon Jackman <acbsjack@gosnell.spc.uchicago.edu>
To: Huiman Xie Barnhart <hxb@panda.sph.emory.edu>
Cc: s-news@utstat.toronto.edu
Subject: Re: Cumulative multivariate normal distribution.
try the toms717 Fortran routine on statlib (under misc, I think).
This was contributed by Alan Miller (a fellow Aussie?!) from CSIRO, I
believe...
header is:
SUBROUTINE MECDF(NDIM, D, RHO, PROB, IER)
C
C This comes from TOMS algorithm 717
C
INTEGER NDIM, IER
DOUBLE PRECISION D(*), PROB, RHO(*)
C-----------------------------------------------------------------
C 6/29/90
C This subroutine is designed to calculate the MVN CDF
C using the Mendell-Elston procedure as described in
C Kamakura (1989). The current version is set up to go
C as high as 19 dimensions (=> 20 MNP alternatives)
C NOTE: Equation (15) in Kamakura has an error.
C
C Specifically, assume that Z is a set of random variables
C with a standard normal distribution with correlations
C stored in RHO (in packed form). Then this subroutine
C calculates Prob[Z(1)>D(1);...; Z(NDIM) > D(NDIM)].
C-----------------------------------------------------------------
I compiled this and put it local.Sqpe (static load only here, ver 3.1;
groan) call it with the following Splus function:
> mvncdf
function(x, rho)
{
# caller function for tom717.f - fortran routine for mecdf
x <- as.vector(x)
ndim <- length(x)
d <- x
# fortran call is
# SUBROUTINE MECDF(NDIM, D, RHO, PROB, IER)
out <- .Fortran("mecdf",
as.integer(ndim),
as.double(d),
as.double(rho),
as.double(1),
as.integer(1))
out
}
regards - simon jackman
========================================================================
Simon Jackman, Dept of Political Science, Univ of Chicago, 5828
S.University Ave, Chicago, IL 60637. ph: 312.702.8075. fax: x.x.1689 "But
the age of chivalry is gone. That of sophisters, economists, and
calculators has succeeded; and the glory of Europe is extinguished for
ever." (Burke, _Reflections on the Revolution in France_) "I can't give
credit enough to the logisticians..." (Norman Schwarzkopf, 2/27/91)
From daniel_heitjan@merck.com
Date: Thu, 02 Mar 1995 07:02 -0500 (EST)
From: Daniel Heitjan <daniel_heitjan@merck.com>
To: Huiman Xie Barnhart <hxb@panda.sph.emory.edu>
Subject: RE: Cumulative multivariate normal distr
There's an Applied Statistics algorithm for the cumulative bivariate normal,
which
can easily be packageded to run in Splus.
Check statlib.
