Wang ZHOU

/*data analysis: 0.95 two sided C.I. by Edgeworth expansions for P(X<Y) for data set (c) in Table 1 of Newcomb (2006b)*/
/*When you apply the following code to your own data set, please do two necessary modifications. One is to change the
sample sizes, "mysizef" and "mysizeg". The 2nd is to input your data sets to "samplef()" and "sampleg()".*/
/*Please note that if all the X are smaller than all the Y, or all the X are bigger than all the Y, the code will not
work since the variances are zero in those two cases. */


#include "math.h"


#define mysizef 5  /* X sample size n1*/
#define mysizeg 5  /* Y sample size n2*/

#define c1 1.644854  /*related to 0.9 nominal*/
#define c2 1.959964  /*related to 0.95 nominal*/

float U2;

float mysamplef[mysizef];
float mysampleg[mysizeg];


float egmysamplef[mysizef]; /*g1(X_i)*/
float egmysampleg[mysizeg]; /*g2(Y_j)*/

void samplef()
{
 mysamplef[0]=1.0; mysamplef[1]=2.0; mysamplef[2]=3.0; mysamplef[3]=4.0;
 mysamplef[4]=6.0;
}

void sampleg()
{
 mysampleg[0]=5.0; mysampleg[1]=7.0; mysampleg[2]=8.0; mysampleg[3]=9.0;
 mysampleg[4]=10.0;
}



float utwo()
{
 int i,j;
 float u=0.0;
 
 for(i=0;i<mysizef;i++)
 {
  for(j=0;j<mysizeg;j++)
  {
   if(mysamplef[i]<mysampleg[j])
   u=u+1.0;
   }
  }
  
  return (u/mysizef/mysizeg);
}



float g1(float x)
{
 int j;
 float u,u1;
 
 u=utwo();
 u1=0.0;

 for(j=0;j<mysizeg;j++)
  {
   if(x<mysampleg[j])
   u1=u1+1.0/mysizeg;
   }
   
 return (u1-u); 
}

float g2(float y)
{
 int i;
 float u,u2;
 
 u=utwo();
 u2=0.0;

 for(i=0;i<mysizef;i++)
  {
   if(mysamplef[i]<y)
   u2=u2+1.0/mysizef;
   }
   
 return (u2-u); 
}
 
float egsample() /*g1(X_1)...g1(X_n1),g2(Y_1)...g2(Y_n2)*/
{
 int i,j;
 float u,v1,v2;
 
 u=utwo();
 for(i=0;i<mysizef;i++)
 {
  egmysamplef[i]=g1(mysamplef[i]);
  }
  
 for(j=0;j<mysizeg;j++)
 {
  egmysampleg[j]=g2(mysampleg[j]);
  }
}


float variance()
{
 int i,j;
 float u,v1,v2;
 
 v1=0.0; v2=0.0;

 
 u=utwo();
 for(i=0;i<mysizef;i++)
 {
   v1=v1+egmysamplef[i]*egmysamplef[i]/(mysizef-1.0)/mysizef;
  }
 
 for(j=0;j<mysizeg;j++)
 {
   v2=v2+egmysampleg[j]*egmysampleg[j]/(mysizeg-1.0)/mysizeg;
  }
  
  return v1+v2;
}

float gamma1()   /*Eg1(X)^3*/
{
 int i,j;
 float u,u1,v1;
 
 v1=0.0;
 
 for(i=0;i<mysizef;i++)
 {
   v1=v1+egmysamplef[i]*egmysamplef[i]*egmysamplef[i]/mysizef;
  }
 
 return v1;
}

float gamma2() /*Eg2(Y)^3*/
{
 int i,j;
 float u,u2,v2;
 
 v2=0.0;

  for(j=0;j<mysizeg;j++)
 {
   v2=v2+egmysampleg[j]*egmysampleg[j]*egmysampleg[j]/mysizeg;
  }  
 
 return v2;
}

float gamma3()
{
 int i,j;
 float v=0.0;
 
 for(i=0;i<mysizef;i++)
 {
  for(j=0;j<mysizeg;j++)
  {
  if(mysamplef[i]<mysampleg[j])
   v=v+egmysamplef[i]*egmysampleg[j];
   }
  }
  
  return (v/mysizef/mysizeg);
}




main()
{
 float eglower2,egupper2;
 float V;
 float hat1,hat2,hat3;
 float d2,N;

 
 N=(mysizef+mysizeg)*1.0;
 
 samplef();
 sampleg();
 egsample();
 
 U2=utwo();
 V=variance();
 printf("U2=%10.6fV=%10.6f\n",U2,V);
 
  hat1=gamma1();
  hat2=gamma2();
  hat3=gamma3();
  d2=-c2*c2*(hat1*N*N/mysizef/mysizef+hat2*N*N/mysizeg/mysizeg)/2.0/sqrt(N)/(N*V)/sqrt(N*V)
    +(c2*c2-1.0)*(hat1*N*N/mysizef/mysizef+hat2*N*N/mysizeg/mysizeg+6.0*hat3*N*N/mysizef/mysizeg)/6.0/sqrt(N)/(N*V)/sqrt(N*V);   

  eglower2=U2-(c2+d2)*sqrt(V);  /*lower bound of CI by Edgeworth*/
  egupper2=U2+(c2-d2)*sqrt(V);  /*upper bound of CI by Edgeworth*/
  
  printf("eglower2=%10.6fegupper2=%10.6f\n",eglower2,egupper2);
}