• Zielgruppen
  • Suche
 

Calculation of orthant normal probabilities for a tridiagonal corelation matrix

Abstract

The following SAS/IML-program calculates orthant probabilities of the multivariate normal distribution possessing a tridiagonal correlation matrix. The  algorithm used here is based on the representations given by Sun (1988) and Sun and Asano (1989). It calculates for a vector containing the elements of the diagonal line adjacent to the main diagonal of the correlation matrix the corresponding multiple integral from -infinity to 0. The generated output shows the values of the sought probability prob. The algorithm works up to dimension 11. As an example, the invokation of the module by
 
 

corr={ 0.7071068 0.5 0.3333333};

prob=orthant(corr);
print prob [format=6.4];

yields the following output:
 

 PROB
0.1364

SAS/IML-Source

start orthant(corr);
  n=ncol(corr)+1;
  r=corr;

  coef=j(1,7,.);
  pai=3.141592653589793238;

  do i=1 to int(n/2)+1;
    coef[i]=1/(2**(n-i+1)*pai**(i-1));
  end;

  if n=1 then pr=coef[1];
  else if n=2  then pr=coef[1]+coef[2]*arsin(r[1]);
  else if n=3  then pr=coef[1]+coef[2]*(arsin(r[1])+arsin(r[2]));
  else if n=4  then pr=coef[1]+coef[2]*sasin(r,4) +coef[3]*fint(r,4);
  else if n=5  then pr=coef[1]+coef[2]*sasin(r,5) +coef[3]*si4(r,5);
  else if n=6  then pr=coef[1]+coef[2]*sasin(r,6) +coef[3]*si4(r,6) +coef[4]*fint(r,6);
  else if n=7  then pr=coef[1]+coef[2]*sasin(r,7) +coef[3]*si4(r,7) +coef[4]*si6(r,7);
  else if n=8  then pr=coef[1]+coef[2]*sasin(r,8) +coef[3]*si4(r,8) +coef[4]*si6(r,8) +coef[5]*fint(r,8);
  else if n=9  then pr=coef[1]+coef[2]*sasin(r,9) +coef[3]*si4(r,9) +coef[4]*si6(r,9) +coef[5]*si8(r,9);
  else if n=10 then pr=coef[1]+coef[2]*sasin(r,10)+coef[3]*si4(r,10)+coef[4]*si6(r,10)+coef[5]*si8(r,10)+coef[6]*fint(r,10);
  else if n=11 then pr=coef[1]+coef[2]*sasin(r,11)+coef[3]*si4(r,11)+coef[4]*si6(r,11)+coef[5]*si8(r,11)+coef[6]*si10(r,11);

  return(pr);
finish;


start sasin(r,n);
  sasin=0;
  do i=1 to n-1;
    sasin=sasin+arsin(r[i]);
  end;
  return(sasin);
finish;


start si4(r,n);
  si4=0;
  r1=dis(r,8,3,n-3);
  do i=1 to n-3;
    p=tran(r1,8,3,i);
    si4=si4+fint(p,4);
  end;
  do i=1 to n-4;
    sum=0;
    do j=i+3 to n-1;
      sum=sum+arsin(r[j]);
    end;
    si4=si4+arsin(r[i])*sum;
  end;
  return(si4);
finish;


start si6(r,n);
  si6=0;
  r1=dis(r,8,3,n-3);
  r2=dis(r,6,5,n-5);
  do i=1 to n-5;
    p=tran(r2,6,5,i);
    si6=si6+fint(p,6);
  end;
  do i=1 to n-6;
    sum=0;
    do j=i+5 to n-1;
      sum=sum+arsin(r[j]);
    end;
    p=tran(r1,8,3,i);
    si6=si6+fint(p,4)*sum;
  end;
  do i=4 to n-3;
    sum=0;
    do j=1 to i-3;
      sum=sum+arsin(r[j]);
    end;
    p=tran(r1,8,3,i);
    si6=si6+fint(p,4)*sum;
  end;
  do i=1 to n-7;
    sum=0;
    do j=i+3 to n-4;
      sum1=0;
      do k=j+3 to n-1;
        sum1=sum1+arsin(r[k]);
      end;
      sum=sum+arsin(r[j])*sum1;
    end;
    si6=si6+arsin(r[i])*sum;
  end;
  return(si6);
finish;


start si8(r,n);
  si8=0;
  r1=dis(r,8,3,n-3);
  r2=dis(r,6,5,n-5);
  r3=dis(r,4,7,n-7);

  do i=1 to n-7;
    p=tran(r3,4,7,i);
    si8=si8+fint(p,8);
  end;
  do i=1 to n-8;
    sum=0;
    do j=i+7 to n-1;
      sum=sum+arsin(r[j]);
    end;
    p=tran(r2,6,5,i);
    si8=si8+fint(p,6)*sum;
  end;
  do i=4 to n-5;
    sum=0;
    do j=1 to i-3;
      sum=sum+arsin(r[j]);
    end;
    p=tran(r2,6,5,i);
    si8=si8+fint(p,6)*sum;
  end;
  do i=1 to n-8;
    sum=0;
    do j=i+5 to n-3;
      p=tran(r1,8,3,j);
      sum=sum+fint(p,4);
    end;
    p=tran(r1,8,3,i);
    si8=si8+fint(p,4)*sum;
  end;
  do i=1 to n-9;
    sum=0;
    do j=i+5 to n-4;
      sum1=0;
      do k=j+3 to n-1;
        sum1=sum1+arsin(r[k]);
      end;
      sum=sum+arsin(r[j])*sum1;
    end;
    p=tran(r1,8,3,i);
    si8=si8+fint(p,4)*sum;
  end;
  do i=7 to n-3;
    sum=0;
    do j=1 to i-6;
      sum1=0;
      do k=j+3 to i-3;
        sum1=sum1+arsin(r[k]);
      end;
      sum=sum+arsin(r[j])*sum1;
    end;
    p=tran(r1,8,3,i);
    si8=si8+fint(p,4)*sum;
  end;
  sum=0;
  do i=9 to n-1;
    sum=sum+arsin(r[i]);
  end;
  p=tran(r1,8,3,4);
  si8=si8+arsin(r[1])*fint(p,4)*sum;
  if n>10 then do;
    sum=0;
    do i=1 to n-9;
      sum=sum+arsin(r[i]);
    end;
    p=tran(r1,8,3,5);
    si8=si8+arsin(r[10])*fint(p,4)*sum;
    si8=si8+arsin(r[1])*arsin(r[4])*arsin(r[7])*arsin(r[10]);
  end;
  return(si8);
finish;


start si10(r,n);
  si10=0;
  r1=dis(r,8,3,n-3);
  r2=dis(r,6,5,n-5);
  r3=dis(r,4,7,n-7);
  r4=dis(r,2,9,n-9);

  do i=1 to n-9;
    p=tran(r4,2,9,i);
    si10=si10+fint(p,10);
  end;
  p=tran(r3,4,7,4);
  q=tran(r3,4,7,1);
  si10=si10+arsin(r[1])*fint(p,8)+arsin(r[10])*fint(q,8);
  p=tran(r1,8,3,1);
  q=tran(r2,6,5,6);
  si10=si10+fint(p,4)*fint(q,6);
  p=tran(r1,8,3,8);
  q=tran(r2,6,5,1);
  si10=si10+fint(p,4)*fint(q,6);
  return(si10);
finish;


start fint(r,m);
  a=j(1,19,.);
  b=j(1,11,.);
  c=j(1,7,.);
  p=j(1,20,.);
  gl=j(2,12,.);

  iint=12;
  g1={-0.981560634246719251 -0.904117256370474857 -0.769902674194304687
      -0.587317954286617447 -0.367831498998180194 -0.125233408511463915,
       0.0471753363865118272 0.106939325995318431  0.160078328543346226
       0.203167426723065922  0.233492536538354809  0.249147045813402785};

  do i=1 to iint/2;
    gl[1,i]=0.5*(g1[1,i]+1);
    gl[1,iint+1-i]=0.5*(-g1[1,i]+1);
    gl[2,i]=0.5*g1[2,i];
    gl[2,iint+1-i]=gl[2,i];
  end;

  i=1;
  do j=1 to m+m-3 by 2;
    p[j]=1;
    p[j+1]=r[i];
    i=i+1;
  end;
  p[m+m-1]=1;

  sum1=0;
  k1=0;
  krit=0;
  do until(k1=iint | krit=1);
    k1=k1+1;
    if m<=8 then do;
      do j=1 to 15;
        a[j]=p[j];
      end;
    end;
    else do;
      ta=1-(p[2]*gl[1,k1])**2;
      a[1]=ta-p[4]*p[4];
      do j=2 to 18 by 2;
        a[j]=p[j+2]*ta;
        a[j+1]=ta;
      end;
    end;

    sum2=0;
    k2=0;
    do until(k2=iint | krit=1);
      k2=k2+1;
      if m<=6 then do;
        do j=1 to 11;
          b[j]=p[j];
        end;
      end;
      else do;
        ub=a[1]*a[3]-(a[2]*gl[1,k2])**2;
        ub1=a[5]*ub;
        b[1]=ub1-a[1]*a[4]*a[4];
        do j=2 to 10 by 2;
          b[j]=a[j+2]*ub;
          b[j+1]=ub1;
        end;
      end;

      sum3=0;
      k3=0;
      do until(k3=iint | krit=1);
        k3=k3+1;
        if m=4 then do;
          do j=1 to 7;
            c[j]=p[j];
          end;
        end;
        else do;
          vc=b[1]*b[3]-(b[2]*gl[1,k3])**2;
          vc1=b[5]*vc;
          c[1]=vc1-b[1]*b[4]*b[4];
          do j=2 to 6 by 2;
            c[j]=b[j+2]*vc;
            c[j+1]=vc1;
          end;
        end;

        sum4=0;
        do k4=1 to iint;
          wd=c[1]*c[3]-(c[2]*gl[1,k4])**2;
          if wd<0 then wd=1E-09;
          wd1=c[5]*wd;
          d11=wd1-c[1]*c[4]*c[4];
          d12=c[6]*wd;
          d22=wd1;
          yyy=d11*d22;
          if yyy<0 then yyy=1E-09;
          xxx=d12/sqrt(yyy);
          if xxx<-1 then xxx=-1+1E-09;
          sum4=sum4+gl[2,k4]/sqrt(wd)*arsin(xxx);
        end;
        sum4=c[2]*sum4;
        if m^=4 then sum3=sum3+gl[2,k3]/sqrt(vc)*sum4;
        else krit=1;
      end;
      if m^=4 then do;
        sum3=b[2]*sum3;
        if m^=6 then sum2=sum2+gl[2,k2]/sqrt(ub)*sum3;
        else krit=1;
      end;
    end;
    if (m^=4 & m^=6) then do;
      sum2=a[2]*sum2;
      if m^=8 then sum1=sum1+gl[2,k1]/sqrt(ta)*sum2;
      else krit=1;
    end;
  end;

  if m=4 then return(sum4);
  else if m=6 then return(sum3);
  else if m=8 then return(sum2);
  else return(p[2]*sum1);
finish;


start dis(r,m,n,k);
  rr=j(k,n,.);
  do i=1 to k;
    do j=1 to n;
      rr[i,j]=r[i-1+j];
    end;
  end;
  return(rr);
finish;


start tran(rr,m,n,k);
  r=j(1,n,.);
  do i=1 to n;
    r[i]=rr[k,i];
  end;
  return(r);
finish;