statistics.h
1 /***************************************************************************
2  statistics.h - description
3  -------------------
4  begin : Fri Apr 6 2001
5  copyright : (C) 2000-2015 by Thies Jochimsen
6  email : thies@jochimsen.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef STATISTICS_H
19 #define STATISTICS_H
20 
21 #include<odindata/data.h>
22 #include<odindata/utils.h>
23 
24 
30 
37 
41  double min;
42 
46  double max;
47 
51  double mean;
52 
56  double stdev;
57 
61  double meandev;
62 
66  friend STD_ostream& operator << (STD_ostream& s, statisticResult stats);
67 
68 };
69 
70 
72 
78 template<typename T, int N_rank>
79 statisticResult statistics(const Array<T,N_rank>& ensemble, const Array<T,N_rank>* mask=0) {
80  Log<OdinData> odinlog("","statistics");
81  statisticResult result;
82 
83  result.min=0.0;
84  result.max=0.0;
85  result.mean=0.0;
86  result.stdev=0.0;
87  result.meandev=0.0;
88 
89 
90  if(mask && !same_shape(ensemble,*mask)) {
91  ODINLOG(odinlog,errorLog) << "size mismatch (ensemble.shape()=" << ensemble.shape() << ") != (mask.shape()=" << mask->shape() << ")" << STD_endl;
92  return result;
93  }
94 
95 
96  int n=ensemble.numElements();
97 
98  Data<T,N_rank> ensemble_copy(ensemble);
99 
100  int nvals=0;
101  TinyVector<int,N_rank> index;
102  for(int i=0; i<n; i++) {
103  index=ensemble_copy.create_index(i);
104  bool include=true;
105  if(mask && (*mask)(index)==0.0) include=false;
106  if(include) {
107  float val=ensemble(index);
108  result.mean+=val;
109  nvals++;
110  if(i==0) {
111  result.min=result.max=val;
112  } else {
113  if(val<result.min) result.min=val;
114  if(val>result.max) result.max=val;
115  }
116  }
117  }
118  result.mean=secureDivision(result.mean,double(nvals));
119 
120  nvals=0;
121  for(int i=0; i<n; i++) {
122  index=ensemble_copy.create_index(i);
123  bool include=true;
124  if(mask && (*mask)(index)==0.0) include=false;
125  if(include) {
126  double diff=result.mean-ensemble(index);
127  result.stdev+=diff*diff;
128  nvals++;
129  }
130  }
131  if(nvals>1) result.stdev=sqrt(result.stdev/double(nvals-1));
132  else result.stdev=0.0;
133 
134  result.meandev=result.stdev/sqrt(double(nvals));
135 
136  return result;
137 }
138 
140 
145 template<typename T, int N_rank>
146 T median(const Array<T,N_rank>& ensemble, const Array<T,N_rank>* mask=0) {
147  T result(0);
148  Data<T,N_rank> ensemble_copy(ensemble);
149  STD_list<T> vallist;
150  int n=ensemble_copy.numElements();
151  if(!n) return result;
152  TinyVector<int,N_rank> index;
153  for(int i=0; i<n; i++) {
154  index=ensemble_copy.create_index(i);
155  bool include=true;
156  if(mask && (*mask)(index)==0.0) include=false;
157  if(include) vallist.push_back(ensemble_copy(index));
158  }
159  vallist.sort();
160 
161  STD_vector<T> vec(list2vector(vallist));
162 
163  if(n%2) {
164  result=vec[(n-1)/2]; // odd number of elements
165  } else {
166  result=0.5*(vec[n/2]+vec[n/2-1]); // even number of elements
167  }
168 
169  return result;
170 }
171 
173 
179 template<typename T, int N_rank>
180 T weightmean(const Array<T,N_rank>& ensemble, const Array<float,N_rank>& weight) {
181  Log<OdinData> odinlog("","weightmean");
182  T result;
183  result=0; // double-stage init for TinyVector
184 
185  Data<T,N_rank> ensemble_copy(ensemble);
186 
187  if(ensemble.shape()!=weight.shape()) {
188  ODINLOG(odinlog,errorLog) << "size mismatch (ensemble.shape()=" << ensemble.shape() << ") != (weight.shape()=" << weight.shape() << ")" << STD_endl;
189  return result;
190  }
191 
192  float weightsum=0.0;
193  TinyVector<int,N_rank> index;
194  for(unsigned int i=0; i<ensemble_copy.numElements(); i++) {
195  index=ensemble_copy.create_index(i);
196  float w=fabs(weight(index));
197  result+=w*ensemble_copy(index);
198  weightsum+=w;
199  }
200  if(weightsum) result/=weightsum;
201 
202  return result;
203 }
204 
205 
209 #endif
210 
TinyVector< int, N_rank > create_index(unsigned long index) const
STD_vector< T > list2vector(const STD_list< T > &src)
Definition: tjvector.h:500
Definition: tjlog.h:218
bool same_shape(const Array< T, N_rank > &a1, const Array< T, N_rank > &a2, const TinyVector< int, N_rank > &dimmask=1)
Definition: utils.h:168
double meandev
Definition: statistics.h:61
friend STD_ostream & operator<<(STD_ostream &s, statisticResult stats)
T weightmean(const Array< T, N_rank > &ensemble, const Array< float, N_rank > &weight)
Definition: statistics.h:180
statisticResult statistics(const Array< T, N_rank > &ensemble, const Array< T, N_rank > *mask=0)
Definition: statistics.h:79
double secureDivision(double numerator, double denominator)
T median(const Array< T, N_rank > &ensemble, const Array< T, N_rank > *mask=0)
Definition: statistics.h:146