40 lines
838 B
C++
40 lines
838 B
C++
#include <cmath>
|
|
#include "nr.h"
|
|
using namespace std;
|
|
|
|
void NR::spear(Vec_I_DP &data1, Vec_I_DP &data2, DP &d, DP &zd, DP &probd,
|
|
DP &rs, DP &probrs)
|
|
{
|
|
int j;
|
|
DP vard,t,sg,sf,fac,en3n,en,df,aved;
|
|
|
|
int n=data1.size();
|
|
Vec_DP wksp1(n),wksp2(n);
|
|
for (j=0;j<n;j++) {
|
|
wksp1[j]=data1[j];
|
|
wksp2[j]=data2[j];
|
|
}
|
|
sort2(wksp1,wksp2);
|
|
crank(wksp1,sf);
|
|
sort2(wksp2,wksp1);
|
|
crank(wksp2,sg);
|
|
d=0.0;
|
|
for (j=0;j<n;j++)
|
|
d += SQR(wksp1[j]-wksp2[j]);
|
|
en=n;
|
|
en3n=en*en*en-en;
|
|
aved=en3n/6.0-(sf+sg)/12.0;
|
|
fac=(1.0-sf/en3n)*(1.0-sg/en3n);
|
|
vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
|
|
zd=(d-aved)/sqrt(vard);
|
|
probd=erfcc(fabs(zd)/1.4142136);
|
|
rs=(1.0-(6.0/en3n)*(d+(sf+sg)/12.0))/sqrt(fac);
|
|
fac=(rs+1.0)*(1.0-rs);
|
|
if (fac > 0.0) {
|
|
t=rs*sqrt((en-2.0)/fac);
|
|
df=en-2.0;
|
|
probrs=betai(0.5*df,0.5,df/(df+t*t));
|
|
} else
|
|
probrs=0.0;
|
|
}
|