Kolmogorov-Smirnov-Test - welche Implementierung für C++
-
Servus,
ich habe ja nun erfahren, dass man oft die GSL nutzt oder nutzen muss, wenn numerische Verfahren gebraucht werden. Wie ist das beim Test einer Datenmenge. Ich brauche eine Implementierung des Kolmogorov-Smirnov-Tests, in R heisst die Funktion ks.test(). Auf den ersten Blick ist der Test nicht Bestandteil der GSL, wo glaubt ihr werde ich fündig ?
Danke vorab, laut Google und Wikipedia fand ich nur eine Quelltext-Funktion für C und eine Umsetzung für C++ - da muss ich doch irgendwas übersehen haben?!
-
Ich denke die Implementierung für C hättest du gefunden? Was brauchst du denn jetzt noch?edit: Oder anders ausgedrückt: Was ist überhaupt deine Frage?
-
Naja, ich dachte es gibt eine copy&paste-Lösungen und mindestens eine verbreitete OpenSource Bibliothek für C++. Ich rechnete mit einer Fülle von gut dokumentierten Ansätzen.
Die C-Implementierung, die ich fand, ist ein PDF aus dem ich den Code abtippen kann, vielleicht fand ich da noch nicht die richtige Adresse, ich meine den Hinweis aus der englischsprachigen Wikipedia:
http://en.wikipedia.org/wiki/Kolmogorov–Smirnov_testWo glaubt ihr werde ich noch fündig?
-
Ich würde da nicht so viel Zeit mit der Suche verbraten. KS-Test sind was, 20 Zeilen? Das hättest du in der Zeit seit du gefragt hast doch längst selber abtippen können. Da ist auch nichts so geheimnisvolles dran, dass du da ein werweißwie durchgestyltes Interface* bräuchtest oder tolle programmiertechnische Tricks anwenden müsstest.
*Ich würde einfach eine Templatefunktion machen die 4 double-Iterator-artige Objekte entgegen nimmt (jeweils für Anfang und Ende) oder 2 Iteratoren und eine Vergleichsfunktion.
-
Hi, ich habe nun mal die C-Variante des Papers auf das in dem englischsprachigen Artikel verwiesen wird hinformatiert und kapiere nun die Bedienung der Funktion nicht. In R gibt man eine Wertemenge und eine Verteilungsform an oder zwei Wertemengen. Hier soll ich einen integer und ein double angeben - wie soll denn das gehen?!
Die Beschaffenheit der Methode und die Erwähnungen im Artikel lassen mich folgern, dass ich da nach externen geprüften Lösungen mich umschauen sollte. Ich muss da lange werkeln und dass ich das fehlerfrei hinbekomme, das ist nicht sicher!
Weiss einer eine solche externe Lösung oder hat eine Methode in der Schublade
?Hier der formatierte Code:
// =============================================================================================================== // The C program for K(n, d) = Pr(Dn < d), aus PDF 'Evaluating Kolmogorov’s Distribution' // The following C program contains the procedure K(n,d), as well as supporting procedures for multiplying // and exponentiating matrices. It is in compact form to save space. To use K(n,d) you need only add a main pro- // gram to a cut-and-paste version of the code listed below. Then make calls to K(n,d) from an int main(){ }. // You should also lead with the usual #include <stdio.h>,#include <math.h> and #include <stdlib.h> void mMultiply(double *A,double *B,double *C,int m) { int i,j,k; double s; for(i=0;i<m;i++) { for(j=0; j<m; j++) { s=0.; for(k=0;k<m;k++) { s+=A[i*m+k]*B[k*m+j]; C[i*m+j]=s; } } } } void mPower(double *A,int eA,double *V,int *eV,int m,int n) { double *B;int eB,i; if(n==1) { for(i=0;i<m*m;i++) { V[i]=A[i];*eV=eA; return; } } mPower(A,eA,V,eV,m,n/2); B=(double*)malloc((m*m)*sizeof(double)); mMultiply(V,V,B,m); eB=2*(*eV); if(n%2==0) { for(i=0;i<m*m;i++) { V[i]=B[i]; *eV=eB; } } else { mMultiply(A,B,V,m); *eV=eA+eB; } if(V[(m/2)*m+(m/2)]>1e140) { for(i=0;i<m*m;i++) { V[i]=V[i]*1e-140;*eV+=140; } } free(B); } double K(int n,double d) { int k,m,i,j,g,eH,eQ; double h,s,*H,*Q; //OMIT NEXT LINE IF YOU REQUIRE >7 DIGIT ACCURACY IN THE RIGHT TAIL s=d*d*n; if(s>7.24||(s>3.76&&n>99)) { return 1-2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s); } k=(int)(n*d)+1; m=2*k-1; h=k-n*d; H=(double*)malloc((m*m)*sizeof(double)); Q=(double*)malloc((m*m)*sizeof(double)); for(i=0;i<m;i++) for(j=0;j<m;j++) { if(i-j+1<0) { H[i*m+j]=0; } else { H[i*m+j]=1; } } for(i=0;i<m;i++) { H[i*m]-=pow(h,i+1); H[(m-1)*m+i]-=pow(h,(m-i)); } H[(m-1)*m]+=(2*h-1>0?pow(2*h-1,m):0); for(i=0;i<m;i++) { for(j=0;j<m;j++) { if(i-j+1>0) { for(g=1;g<=i-j+1;g++) { H[i*m+j]/=g; } } } } eH=0; mPower(H,eH,Q,&eQ,m,n); s=Q[(k-1)*m+k-1]; for(i=1;i<=n;i++) { s=s*i/n; if(s<1e-140) { s*=1e140; eQ-=140; } } s*=pow(10.,eQ); free(H); free(Q); return s; } // ==== Ende KS-Test ===========================================================
-
Jay1980 schrieb:
Weiss einer eine solche externe Lösung oder hat eine Methode in der Schublade
?Numerical Recipes enthält dafür einen kurzen, kompakten Code mit offensichtlicher Art und Weise wie er zu bedienen ist.
-
Erstmal danke, ich habe gerade in der dritten Auflage des Buches nachgesehen, ja da ist es möglich nachzulesen, was da passiert (S. 737 ff). Leider sind aber dort nur die Header-Dateien sichtbar und diese greifen auf Objekte von andernorts zu.
Dazu kommt die Lizenzsache, ich habe gelesen, dass NR-Code recht schwer einsetzbar ist. Später muss ich zwar nur einmal was 'für mich' berechnen. Ich schau nun, wie komplex die eingesetzten Objekte sind, aber zusammenfassend ist der NR-Code-Snippet weit weg für mich vom 'kurz-kompakt-verständlich-aufseigeneproblemumlegbar' ... mal schauen.
Edit: Okay ... der Aufbau der verwendeten Objekte wird andernorts im Buch beschrieben, ggf. geht da doch was!
-
Danke, die ersten kleinen NR-Beispiele laufen, ich denke ich werde noch oft auf dieses Buch zurückgreifen - zwar etwas Gefrickel aber danach weiss man wenigstens was da wie funktioniert.
-
Ich habe doch noch Probleme mit der Numerical-Recipes-Variante, die Werte sind komisch. Test 1 und Test 2 (die Tests stehen ganz unten oberhalb der main-Methode) liefern mir p-Werte um die 0.99 also sehr hoch, die gleichen Werte, wenn diese Werte viele Nachkommastellen haben, liefert mir sehr kleine Werte. Beim KS-Test geht es ja um die größte Distanz zwischen den Wertepaaren, diese ist beim dritten Test einfach sehr sehr hoch, bei 0.7 ...
Wie man sieht übergebe ich die Werte die ich aus einer kumulierten Funktion, ich dachte der KS-Test von NR will das so (S. 738ff Numerical Recipes, dritte Edition). Danke vorab, falls einer die NR-Lib zur Hand hat und mal nachtesten will oder ein kleines Beispiel parat hat, dass die KS-Test Anwendung bei ihm zeigt.
Hier ist das Code-Ausschnitt, er ist natürlich nur kompilierbar, wenn man die header-Dateien verfügbar hat, also die NR-Bibliothek.
#include "nr3.h" #include "sort.h" // NR-Struktur fuer KS-Verteilung struct KSdist { // Kap 14.3.3. S.737 + Kap 6.14.56 S. 334 double pks(double z) { // Return cumulative distribution function. if (z < 0.) throw ("bad z in KSdist"); if (z == 0.) return 0.; if (z < 1.18) { double y = exp(-1.23370055013616983 / SQR(z)); return 2.25675833419102515 * sqrt(-log(y))*(y + pow(y, 9) + pow(y, 25) + pow(y, 49)); } else { double x = exp(-2. * SQR(z)); return 1. - 2. * (x - pow(x, 4) + pow(x, 9)); } } // Return complementary cumulative distribution function. double qks(double z) { if (z < 0.) throw ("bad z in KSdist"); if (z == 0.) return 1.; if (z < 1.18) return 1. - pks(z); double x = exp(-2. * (z * z)); return 2. * (x - pow(x, 4) + pow(x, 9)); } // more code in the book ... }; // NR-Funkton zum KS-Test void kstwo(VecDoub_IO &data1, VecDoub_IO &data2, Doub &d, Doub &prob) { // Given an array data1[0..n1-1], and an array data2[0..n2-1], this routine returns the K–S // statistic d and the p-value prob for the null hypothesis that the data sets are drawn from the // same distribution. Small values of prob show that the cumulative distribution function of data1 // is significantly different from that of data2. The arrays data1 and data2 are modified by being // sorted into ascending order. int j1 = 0, j2 = 0, n1 = data1.size(), n2 = data2.size(); double d1, d2, dt, en1, en2, en, fn1 = 0.0, fn2 = 0.0; KSdist ks; sort(data1); sort(data1); // Ausgabe der beiden sortierten Collections //int nrRealVecLength = data1.size(); //int nrGumbelVecLength = data2.size(); //for (int i = 0; i < nrRealVecLength; i++) { // cout << "Debug nrRealVec data1 ... Lauf " << i << ": " << data1[i] << endl; //} //for (int i = 0; i < nrGumbelVecLength; i++) { // cout << "Debug nrGumbelVec data2 ... Lauf " << i << ": " << data2[i] << endl; //} en1 = n1; en2 = n2; d = 0.0; while (j1 < n1 && j2 < n2) // If we are not done... { if ((d1 = data1[j1]) <= (d2 = data2[j2])) // Next step is in data1. { do { fn1 = ++j1 / en1; } while (j1 < n1 && d1 == data1[j1]); } if (d2 <= d1) // Next step is in data2. { do { fn2 = ++j2 / en2; } while (j2 < n2 && d2 == data2[j2]); } if ((dt = abs(fn2 - fn1)) > d) { //cout << " fn2 " << fn2 << "; fn1 " << fn1 << "; d " << d << endl; d = dt; } } en = sqrt(en1 * en2 / (en1 + en2)); //cout << "* en " << en << " aus en1 " << en1 << " und en2 " << en2 << endl; prob = ks.qks( (en + 0.12 + 0.11 / en) * d); // Compute p-value. } void testkstwo(){ // Test 1 double nrDis = 0.0; double nrProb = 0.0; int elementsToCompare = 100; // nutzlos, da bei Direktinitialisierung die Größe zur Kompilierzeit da sein muss VecDoub vOne(100); VecDoub vTwo(100); double v1[100] = { 0, 0, 0, 0, 0, 24, 98, 196, 329, 476, 610, 730, 819, 888, 934, 965, 981, 988, 995, 998, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; double v2[100] = { 0, 0, 0, 0, 0, 4, 31, 111, 250, 418, 577, 707, 804, 872, 917, 947, 966, 979, 986, 991, 995, 997, 998, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; for ( int i = 0; i < elementsToCompare; i++ ) { vOne[i] = v1[i]; vTwo[i] = v2[i]; cout << " -> Lauf - " << i << "; Vergleich an vOne - " << vOne[i] << "; vTwo " << vTwo[i] << "; " << endl; } kstwo( vOne, vTwo, nrDis, nrProb ); cout << "KS-Test: Dist: " << nrDis << "; P: " << nrProb << endl; // Test 2 double nrDisT2 = 0.0; double nrProbT2 = 0.0; int elementsToCompareT2 = 100; // nutzlos, da bei Direktinitialisierung die Größe zur Kompilierzeit da sein muss VecDoub vOneT2(100); VecDoub vTwoT2(100); double v1T2[100] = { 0, 0, 0, 0, 0, 24, 98, 196, 329, 476, 610, 730, 819, 888, 934, 965, 981, 988, 995, 998, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; double v2T2[100] = { 0, 0, 0, 0, 0, 4, 31, 111, 250, 418, 577, 707, 804, 872, 917, 947.56, 966, 979.88, 986, 991, 995, 997, 998, 999, 999, 999, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; for ( int i = 0; i < elementsToCompare; i++ ) { vOneT2[i] = v1T2[i]; vTwoT2[i] = v2T2[i]; cout << " -> Lauf - " << i << "; Vergleich an vOneT2 - " << vOneT2[i] << "; vTwoT2 " << vTwoT2[i] << "; " << endl; } kstwo( vOneT2, vTwoT2, nrDisT2, nrProbT2 ); cout << "KS-Test: DistT2: " << nrDisT2 << "; PT2: " << nrProbT2 << endl; // Test 3 double nrDisT3 = 0.0; double nrProbT3 = 0.0; int elementsToCompareT3 = 100; // nutzlos, da bei Direktinitialisierung die Größe zur Kompilierzeit da sein muss VecDoub vOneT3(100); VecDoub vTwoT3(100); double v1T3[100] = { 0, 0, 0, 0, 0, 24, 98, 196, 329, 476, 610, 730, 819, 888, 934, 965, 981, 988, 995, 998, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; double v2T3[100] = { 6.24211e-22, 5.61377e-13, 2.45976e-07, 0.000882158, 0.153194, 3.94838, 30.589, 111.12, 250.483, 418.009, 577.19, 707.315, 803.979, 871.56, 917.028, 946.887, 966.198, 978.567, 986.442, 991.436, 994.595, 996.591, 997.851, 998.645, 999.146, 999.462, 999.661, 999.786, 999.865, 999.915, 999.947, 999.966, 999.979, 999.987, 999.992, 999.995, 999.997, 999.998, 999.999, 999.999, 999.999, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000 }; double flooredV2T3[100]; for ( int i = 0; i < elementsToCompare; i++ ) { flooredV2T3[i] = floor( v2T3[i] ); cout << "Rundungsalarm - Lauf " << i << " : " << v2T3[i] << " zu " << flooredV2T3[i] << endl; } for ( int i = 0; i < elementsToCompare; i++ ) { vOneT3[i] = v1T3[i]; vTwoT3[i] = flooredV2T3[i]; cout << " -> Lauf - " << i << "; Vergleich an vOneT3 - " << vOneT3[i] << "; vTwoT3 " << vTwoT3[i] << "; " << endl; } kstwo( vOneT3, vTwoT3, nrDisT3, nrProbT3 ); cout << "KS-Test: DistT3: " << nrDisT3 << "; PT3: " << nrProbT3 << endl; } int main(int argc, char** argv) { testkstwo(); return 0; }
-
Falls jemand den Test kennt und dazu ein kleines Zahlenbeispiel posten kann, also seine Werte und sein Ergebnis, dann kann ich das mit den von mir kommenden Ergebnissen vergleichen.
