aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-07-12 10:06:47 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:53 +0100
commit6264ff6e0ac13fcc22afcb4ce9f98a920444a92c (patch)
tree4a8d7fd48fbb1b64815da59faf29e921a73d7542
parentf672ed3da7d82b5519b1b5aa77c82e16cc11e0a5 (diff)
compare_hkl: Add Pearson correlation
-rw-r--r--src/compare_hkl.c4
-rw-r--r--src/statistics.c36
-rw-r--r--src/statistics.h3
3 files changed, 42 insertions, 1 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 3e313733..0dffd3a0 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
char *afile = NULL;
char *bfile = NULL;
signed int h, k, l;
- double scale, R2, Rmerge;
+ double scale, R2, Rmerge, pearson;
unsigned int *c1;
unsigned int *c2;
int i;
@@ -145,6 +145,8 @@ int main(int argc, char *argv[])
STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
+ pearson = stat_pearson(ref1, c1, ref2, c2);
+ STATUS("Pearson r = %5.4f\n", pearson);
if ( outfile != NULL ) {
write_reflections(outfile, NULL, out, NULL, 1, cell, 1);
diff --git a/src/statistics.c b/src/statistics.c
index da3d95c4..f3fb21c0 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -18,6 +18,7 @@
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>
+#include <gsl/gsl_statistics.h>
#include "statistics.h"
#include "utils.h"
@@ -249,3 +250,38 @@ double stat_r2(const double *ref1, const unsigned int *c1,
{
return r_minimised(ref1, c1, ref2, c2, scalep, R_2);
}
+
+
+double stat_pearson(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2)
+{
+ double vec1[4096];
+ double vec2[4096];
+ signed int h, k, l;
+ int i = 0;
+
+ for ( l=-INDMAX; l<INDMAX; l++ ) {
+ for ( k=-INDMAX; k<INDMAX; k++ ) {
+ for ( h=-INDMAX; h<INDMAX; h++ ) {
+
+ double i1, i2;
+ unsigned int c1s, c2s;
+
+ c1s = lookup_count(c1, h, k, l);
+ c2s = lookup_count(c2, h, k, l);
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ i2 = lookup_intensity(ref2, h, k, l);
+
+ if ( c1s && c2s && (i1>0.0) && (i2>0.0) ) {
+ vec1[i] = sqrt(i1 / (double)c1s);
+ vec2[i] = sqrt(i2 / (double)c2s);
+ i++;
+ }
+
+ }
+ }
+ }
+
+ return gsl_stats_correlation(vec1, 1, vec2, 1, i);
+}
diff --git a/src/statistics.h b/src/statistics.h
index e19d4b29..5e3ca52c 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -27,4 +27,7 @@ extern double stat_rmerge(const double *ref1, const unsigned int *c1,
const double *ref2, const unsigned int *c2,
double *scalep);
+extern double stat_pearson(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2);
+
#endif /* STATISTICS_H */