作者:

标题:

C++ code:

// A simulation sample
//
// The returned value "sum" calculates the test statistics Dnq before standardization.
// The returned value "Sigma" calculates the variance of the test statistics Dnq.
#include <iostream>
#include <fstream>
#include <Eigen/Dense>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
using namespace Eigen;
using namespace std;
#define PI 3.141592654

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
/*Uniform0-1, idum negative*/
double ran1(long *idum)
{
int j;
long k;
static long iy=0;
static long iv[NTAB];
double temp;
if (*idum <= 0 || !iy) {  
    if (-(*idum) < 1) *idum=1; 
    else *idum = -(*idum);
    for (j=NTAB+7;j>=0;j--) { 
        k=(*idum)/IQ;
        *idum=IA*(*idum-k*IQ)-IR*k;
        if (*idum < 0) *idum += IM;
        if (j < NTAB) iv[j] = *idum;
    }
    iy=iv[0];
}
k=(*idum)/IQ; 
*idum=IA*(*idum-k*IQ)-IR*k;
if (*idum < 0) *idum += IM;
j=iy/NDIV; 
iy=iv[j]; 
iv[j] = *idum; 
if ((temp=AM*iy) > RNMX) return RNMX; 
else return temp;
}
/*Uniform0-1, idum negative*/
/*N(0,1)*/
double gasdev(long *idum)
{
 double ran1(long *idum);
 static int iset=0;
 static double gset;
 double fac,rsq,v1,v2;
 if (*idum<0) iset=0;
 if(iset == 0){
  do{
   v1=2.0*ran1(idum)-1.0;
   v2=2.0*ran1(idum)-1.0;
   rsq=v1*v1+v2*v2;
  }while(rsq >= 1.0 || rsq == 0.0);
  fac=sqrt(-2.0*log(rsq)/rsq);
  gset=v1*fac;
  iset=1;
  return v2*fac;
 }else{
  iset=0;
  return gset;
 }
}
/*N(0,1)*/

/*Gamma(int ia,1),ia>=1.*/
double gamdev(int ia,long *idum)
{
 double ran1(long *idum);
 int j;
 double am,e,s,v1,v2,x,y;
 if (ia<6) {
  x=1.0;
  for (j=1;j<=ia;j++) x*=ran1(idum);
  x=-log(x);
 }else {
  do{
   do{
    do{
     v1=ran1(idum);
     v2=2.0*ran1(idum)-1.0;
    }while (v1*v1+v2*v2>1.0);
    y=v2/v1;
    am=ia-1;
    s=sqrt(2.0*am+1.0);
    x=s*y+am;
   }while (x<=0.0);
   e=(1.0+y*y)*exp(am*log(x/am)-s*y);
  }while (ran1(idum)>e);
 }
 return x;
}
/*Gamma(int ia,1),ia>=1.*/

double P(int i,int j)
{
 double ans=1;
 int i1;
 for(i1=1;i1<=j;i1++)
  ans=ans*i1;
 for(i1=1;i1<=j-i;i1++)
  ans=ans/double(i1);
 return ans;
}

int main()
{
 int n,p,i,j,l,rep,v,k0,q,qnum; /*n: sample size; p: dimension; k0: true bandwidth of Sigma */
 n=50;
 p=50;
 rep=1000;
 k0=5;
 qnum=49;
 MatrixXd G,z(n,p),x(n,p),sum(rep,(qnum+1)),sigma(rep,(qnum+1)),O,C,xcolsum,xmean;
 G = MatrixXd::Identity(p,p);
 O = MatrixXd::Ones(n,n);
 for(j=1;j<(k0+1);j++)
 {
  for(i=0;i<p-j;i++)
  {
   G(i+j,i) = 0.4;
  }
 }/* Generate Gamma matrix;*/
 /*calculate the excution time*/
 clock_t tStart = clock();
 for (v=0;v<rep;v++)
 {
  long idum=-(v+12345);
  for (i=0;i<n;i++)
  {
   for (j=0;j<p;j++)
   {
    z(i,j)=gasdev(&idum);
   }
  }
  
  x=z*G.transpose();
  xcolsum = x.colwise().sum();
  xmean = O*x/n;
  C=(x-O*x/n).transpose()*(x-O*x/n)/(n-1);
  /*Sample Covariance Matrix of X*/

  for (q=0;q<(qnum+1);q++)
  {
   MatrixXd Y(n,(p-q));
   for (i=0;i<n;i++)
   {
    for (l=0;l<(p-q);l++)
    {
     Y(i,l) = (x(i,l)-xmean(l))*(x(i,(l+q))-xmean(l+q)) - C(l,(l+q));
    }
   }
   MatrixXd YY(n,n),YYcolsum;
   YY=Y*Y.transpose();

   double sigma1,sum1,temp1,temp2,temp3;
   sigma1 = 0;
   sum1=0;
   temp1=0;
   temp2=0;
   temp3=0;
   for(i=0;i<n;i++)
   {
    for(j=0;j<i;j++)
    {
     for (l=0;l<p-q;l++)
     {
      temp1 += x(i,l)*x(i,l+q)*x(j,l)*x(j,l+q);
      temp2 += (x(i,l)*x(j,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(j,l)*x(j,l+q)*xcolsum(l+q));
      temp3 += (x(i,l)*x(j,l+q)*xcolsum(l)*xcolsum(l+q) + x(i,l)*x(i,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(i,l)*x(j,l+q)*xcolsum(l+q) - x(i,l)*x(j,l+q)*x(j,l+q)*xcolsum(l));
     }
    } 
    for(j=(i+1);j<n;j++)
    {
     for (l=0;l<p-q;l++)
     {
      temp1 += x(i,l)*x(i,l+q)*x(j,l)*x(j,l+q);
      temp2 += (x(i,l)*x(j,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(j,l)*x(j,l+q)*xcolsum(l+q));
      temp3 += (x(i,l)*x(j,l+q)*xcolsum(l)*xcolsum(l+q) + x(i,l)*x(i,l)*x(j,l+q)*x(j,l+q) - x(i,l)*x(i,l)*x(j,l+q)*xcolsum(l+q) - x(i,l)*x(j,l+q)*x(j,l+q)*xcolsum(l));
     }
    } 
    sigma1 += YY(i,i)*YY(i,i);
   } 
   sum(v,q)=temp1/n/(n-3)+temp2*(2*n-3)/n/(n-1)/(n-2)/(n-3)+temp3/n/(n-1)/(n-2)/(n-3);
   sigma(v,q)=((YY*YY).trace()-sigma1)/n/(n-1);
  }
 }

 ofstream f;
 f.open ("Sigma_n50p50normal.txt");
 f << sigma;
 f.close();
 ofstream f1;
 f1.open ("Dnq_n50p50normal.txt");
 f1 << sum;
 f1.close();
 return 0;
}