Add computation of the AD for different freq ranges

parent 6c804ffa
......@@ -3339,7 +3339,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
FILE *pFile;
FILE *cFile;
int N,n,nn;
int N,n,nn, i;
double f,dre,dim,Sn,Sn25,Sn75,Sn05,Sn95;
int Nsamp;
......@@ -3508,38 +3508,85 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
double S,A;
// Anderson-Darling statistic
S=0.0;
double lx,ly;
//We want to compute AD for different frequency ranges, from srate 256 up to the actual srate of the run.
for(n=0; n<Nsamp; n++)
double factorAD = ((Nsamp+2.*fmin*Tobs)/256./Tobs);// srate/256
int iAD = (int)(log(factorAD)/log(2.));
double p, c;
int newNsamp;
//We also need the unsorted samples again
sprintf(filename,"%s/fourier_domain_%s_median_residual_%s.dat",outdir,model,ifos[ifo]);
dFile = fopen(filename,"r");
sprintf(filename,"%s/%s_median_PSD_%s.dat",outdir,model,ifos[ifo]);
nFile = fopen(filename,"r");
nn=0;
for(n=0; n<N; n++)
{
double x = gsl_cdf_ugaussian_P(samples[n]);
if(x<0.999)
{
lx = log(x);
ly = log(1.-x);
}
else
fscanf(dFile,"%lg %lg %lg",&f,&dre,&dim);
fscanf(nFile,"%lg %lg %lg %lg %lg %lg",&f,&Sn,&Sn25,&Sn75,&Sn05,&Sn95);
if(f>=fmin)
{
double u = samples[n];
double z = 1.0-1.0/(u*u)+3.0/(u*u*u*u)-15.0/(u*u*u*u*u*u)+105.0/(u*u*u*u*u*u*u*u);
double y = z*exp(-u*u/2.0)/(u*sqrt(M_PI));
lx = -y;
ly = -u*u/2.0 +log(z/(u*sqrt(M_PI)));
//samples[2*nn] = dre/sqrt(Sn/(Tobs/2.));
//samples[2*nn+1] = dim/sqrt(Sn/(Tobs/2.));
samples[2*nn] = dre/sqrt(Sn/2.);//the factor of 2 works, but we need to understand why better
samples[2*nn+1] = dim/sqrt(Sn/2.);
nn++;
}
S += ((2.*(double)(n+1)-1.0)*lx + ((double)(2*Nsamp)-2.*(double)(n+1)+1.0)*ly)/(double)(Nsamp);
}
A = (double)(-1*Nsamp) - S;
// p-value of A^2
double p = 1.-AD(Nsamp,A);
free(samples);
fclose(dFile);
fclose(nFile);
sprintf(filename,"%s/%s_anderson_darling_%s.dat",outdir,model,ifos[ifo]);
FILE *afile = fopen(filename,"w");
fprintf(afile,"%lg %lg\n",A, p);
c=1;
for (i=0; i<=iAD; i++){
newNsamp = (int)(Nsamp+2.*fmin*Tobs)/(c)-2.*fmin*Tobs;// Samples only up to srate = (Nsamp+2.*fmin*Tobs)/(2**i)/Tobs, namely 256, 512, ...
S=0.0;
double *newsamples=NULL;
newsamples=malloc(newNsamp*sizeof(double));
for(n=0; n<newNsamp; n++) newsamples[n]=samples[n];
gsl_sort(newsamples,1,newNsamp);
for(n=0; n<newNsamp; n++)
{
double x = gsl_cdf_ugaussian_P(newsamples[n]);
if(x<0.999)
{
lx = log(x);
ly = log(1.-x);
}
else
{
double u = newsamples[n];
double z = 1.0-1.0/(u*u)+3.0/(u*u*u*u)-15.0/(u*u*u*u*u*u)+105.0/(u*u*u*u*u*u*u*u);
double y = z*exp(-u*u/2.0)/(u*sqrt(M_PI));
lx = -y;
ly = -u*u/2.0 +log(z/(u*sqrt(M_PI)));
}
S += ((2.*(double)(n+1)-1.0)*lx + ((double)(2*newNsamp)-2.*(double)(n+1)+1.0)*ly)/(double)(newNsamp);
}
A = (double)(-1*newNsamp) - S;
// p-value of A^2
p = 1.-AD(newNsamp,A);
fprintf(afile,"%g %lg %lg\n",(Nsamp+2.*fmin*Tobs)/(c)/Tobs, A, p);
c *=2;
free(newsamples);
}
free(samples);
fclose(afile);
// Besides the median PSD, we also compute the AD and p value with the burnin PSD
double *samplesB=NULL;
......@@ -3575,6 +3622,7 @@ void anderson_darling_test(char *outdir, int ifo, double fmin, double Tobs, char
gsl_sort(samplesB,1,Nsamp);
for(n=0; n<Nsamp; n++)
{
double x = gsl_cdf_ugaussian_P(samplesB[n]);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment