diff options
| author | Thomas White <taw@physics.org> | 2010-07-12 10:06:47 +0200 |
|---|---|---|
| committer | Thomas White <taw@physics.org> | 2012-02-22 15:26:53 +0100 |
| commit | 6264ff6e0ac13fcc22afcb4ce9f98a920444a92c (patch) | |
| tree | 4a8d7fd48fbb1b64815da59faf29e921a73d7542 | |
| parent | f672ed3da7d82b5519b1b5aa77c82e16cc11e0a5 (diff) | |
compare_hkl: Add Pearson correlation
| -rw-r--r-- | src/compare_hkl.c | 4 | ||||
| -rw-r--r-- | src/statistics.c | 36 | ||||
| -rw-r--r-- | src/statistics.h | 3 |
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 */ |