SNAP Library 4.0, User Reference  2017-07-27 13:18:06
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
xmath.cpp
Go to the documentation of this file.
1 // Mathematical-Utilities
3 double TMath::E=2.71828182845904523536;
4 double TMath::Pi=3.14159265358979323846;
5 double TMath::LogOf2=log(double(2));
6 
8 // Special-Functions
10  double& gamser, const double& a, const double& x, double& gln){
11  static const int ITMAX=100;
12  static const double EPS=3.0e-7;
13  int n;
14  double sum, del, ap;
15 
16  gln=LnGamma(a);
17  if (x <= 0.0){
18  IAssert(x>=0); /*if (x < 0.0) nrerror("x less than 0 in routine gser");*/
19  gamser=0.0;
20  return;
21  } else {
22  ap=a;
23  del=sum=1.0/a;
24  for (n=1; n<=ITMAX; n++){
25  ++ap;
26  del *= x/ap;
27  sum += del;
28  if (fabs(del) < fabs(sum)*EPS){
29  gamser=sum*exp(-x+a*log(x)-(gln));
30  return;
31  }
32  }
33  Fail; /*nrerror("a too large, ITMAX too small in routine gser");*/
34  return;
35  }
36 }
37 
39  double& gammcf, const double& a, const double& x, double& gln){
40  static const int ITMAX=100;
41  static const double EPS=3.0e-7;
42  static const double FPMIN=1.0e-30;
43  int i;
44  double an, b, c, d, del, h;
45 
46  gln=LnGamma(a);
47  b=x+1.0-a;
48  c=1.0/FPMIN;
49  d=1.0/b;
50  h=d;
51  for (i=1;i<=ITMAX;i++){
52  an = -i*(i-a);
53  b += 2.0;
54  d=an*d+b;
55  if (fabs(d) < FPMIN) d=FPMIN;
56  c=b+an/c;
57  if (fabs(c) < FPMIN) c=FPMIN;
58  d=1.0/d;
59  del=d*c;
60  h *= del;
61  if (fabs(del-1.0) < EPS) break;
62  }
63  IAssert(i<=ITMAX);
64  /*if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");*/
65  gammcf=exp(-x+a*log(x)-(gln))*h;
66 }
67 
68 double TSpecFunc::GammaQ/*gammq*/(const double& a, const double& x){
69  IAssert((x>=0)&&(a>0));
70  double gamser, gammcf, gln;
71  if (x<(a+1.0)){
72  GammaPSeries(gamser,a,x,gln);
73  return 1.0-gamser;
74  } else {
75  GammaQContFrac(gammcf,a,x,gln);
76  return gammcf;
77  }
78 }
79 
80 double TSpecFunc::LnGamma/*gammln*/(const double& xx){
81  double x, y, tmp, ser;
82  static double cof[6]={76.18009172947146,-86.50532032941677,
83  24.01409824083091,-1.231739572450155,
84  0.1208650973866179e-2,-0.5395239384953e-5};
85  int j;
86 
87  y=x=xx;
88  tmp=x+5.5;
89  tmp -= (x+0.5)*log(tmp);
90  ser=1.000000000190015;
91  for (j=0;j<=5;j++) ser += cof[j]/++y;
92  return -tmp+log(2.5066282746310005*ser/x);
93 }
94 
95 double TSpecFunc::LnComb(const int& n, const int& k){
96  return LnGamma(n+1)-LnGamma(k+1)-LnGamma(n-k+1);
97 }
98 
99 double TSpecFunc::BetaCf(const double& a, const double& b, const double& x){
100  static const double MAXIT=100;
101  static const double EPS=3.0e-7;
102  static const double FPMIN=1.0e-30;
103  int m,m2;
104  double aa,c,d,del,h,qab,qam,qap;
105 
106  qab=a+b;
107  qap=a+1.0;
108  qam=a-1.0;
109  c=1.0;
110  d=1.0-qab*x/qap;
111  if (fabs(d) < FPMIN) d=FPMIN;
112  d=1.0/d;
113  h=d;
114  for (m=1;m<=MAXIT;m++) {
115  m2=2*m;
116  aa=m*(b-m)*x/((qam+m2)*(a+m2));
117  d=1.0+aa*d;
118  if (fabs(d) < FPMIN) d=FPMIN;
119  c=1.0+aa/c;
120  if (fabs(c) < FPMIN) c=FPMIN;
121  d=1.0/d;
122  h *= d*c;
123  aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
124  d=1.0+aa*d;
125  if (fabs(d) < FPMIN) d=FPMIN;
126  c=1.0+aa/c;
127  if (fabs(c) < FPMIN) c=FPMIN;
128  d=1.0/d;
129  del=d*c;
130  h *= del;
131  if (fabs(del-1.0) < EPS) break;
132  }
133  if (m > MAXIT){Fail;}// a or b too big, or MAXIT too small in betacf
134  return h;
135 }
136 
137 double TSpecFunc::BetaI(const double& a, const double& b, const double& x){
138  double bt;
139 
140  if (x < 0.0 || x > 1.0){Fail;} // Bad x in routine betai
141  if (x == 0.0 || x == 1.0) bt=0.0;
142  else
143  bt=exp(LnGamma(a+b)-LnGamma(a)-LnGamma(b)+a*log(x)+b*log(1.0-x));
144  if (x < (a+1.0)/(a+b+2.0))
145  return bt*BetaCf(a,b,x)/a;
146  else
147  return 1.0-bt*BetaCf(b,a,1.0-x)/b;
148 }
149 
151  const TVec<TFltPr>& XY, double& A, double& B,
152  double& SigA, double& SigB, double& Chi2, double& R2) {
153  // y = a + bx :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
154  int i;
155  double t, sxoss, sx = 0.0, sy = 0.0, st2 = 0.0, ss, sigdat;
156 
157  A = B = SigA = SigB = Chi2 = 0.0;
158  for (i = 0; i < XY.Len(); i++) {
159  sx += XY[i].Val1;
160  sy += XY[i].Val2;
161  }
162  ss = XY.Len();
163  sxoss = sx / ss;
164  for (i = 0; i <XY.Len(); i++) {
165  t = XY[i].Val1 - sxoss;
166  st2 += t*t;
167  B += t * XY[i].Val2;
168  }
169  B /= st2;
170  A = (sy - sx * B) / ss;
171  SigA = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
172  SigB = sqrt(1.0 / st2);
173  for (i = 0; i < XY.Len(); i++)
174  Chi2 += TMath::Sqr(XY[i].Val2 - A - B * XY[i].Val1);
175  sigdat = sqrt(Chi2 / (XY.Len() - 2));
176  SigA *= sigdat;
177  SigB *= sigdat;
178 
179  // calculate R squared
180  { double N = XY.Len(), sXY=0.0, sX=0.0, sY=0.0, sSqX=0.0, sSqY=0.0;
181  for (int s =0; s < XY.Len(); s++) {
182  sX += XY[s].Val1; sY += XY[s].Val2;
183  sXY += XY[s].Val1 * XY[s].Val2;
184  sSqX += TMath::Sqr(XY[s].Val1);
185  sSqY += TMath::Sqr(XY[s].Val2);
186  }
187  R2 = TMath::Sqr(N*sXY - sX*sY) / ((N*sSqX - sX*sX) * (N*sSqY - sY*sY)); }
188  if (1.1 < R2 || -1.1 > R2) R2 = 0.0;
189  if (_isnan(A) || ! _finite(A)) A = 0.0;
190  if (_isnan(B) || ! _finite(B)) B = 0.0;
191 }
192 
193 void TSpecFunc::PowerFit(const TVec<TFltPr>& XY, double& A, double& B,
194  double& SigA, double& SigB, double& Chi2, double& R2) {
195  // y = a x^b :: SigA (SigB): A's uncertainty; Chi2: std dev on all points
196  // log fit
197  double AA, BB;
198  TFltPrV LogXY(XY.Len(), 0);
199  for (int s = 0; s < XY.Len(); s++) {
200  LogXY.Add(TFltPr(log((double)XY[s].Val1), log((double)XY[s].Val2)));
201  }
202  TSpecFunc::LinearFit(LogXY, AA, BB, SigA, SigB, Chi2, R2);
203  A = exp(AA); B = BB;
204  if (_isnan(AA) || ! _finite(AA)) A = 0.0;
205  if (_isnan(BB) || ! _finite(BB)) B = 0.0;
206 }
207 
208 void TSpecFunc::LogFit(const TVec<TFltPr>& XY, double& A, double& B,
209  double& SigA, double& SigB, double& Chi2, double& R2) {
210  // Y = A + B*log(X)
211  TFltPrV LogXY(XY.Len(), 0);
212  for (int s = 0; s < XY.Len(); s++) {
213  LogXY.Add(TFltPr(log((double)XY[s].Val1), XY[s].Val2));
214  }
215  TSpecFunc::LinearFit(LogXY, A, B, SigA, SigB, Chi2, R2);
216 }
217 
218 void TSpecFunc::ExpFit(const TVec<TFltPr>& XY, double& A, double& B,
219  double& SigA, double& SigB, double& Chi2, double& R2) {
220  // Y = A * exp(B*X)
221  TFltPrV XLogY(XY.Len(), 0);
222  double AA, BB;
223  for (int s = 0; s < XY.Len(); s++) {
224  XLogY.Add(TFltPr(XY[s].Val1, log((double)XY[s].Val2)));
225  }
226  TSpecFunc::LinearFit(XLogY, AA, BB, SigA, SigB, Chi2, R2);
227  A = exp(AA);
228  B = BB;
229 }
230 
231 double TSpecFunc::Entropy(const TIntV& ValV) {
232  TFltV NewValV(ValV.Len());
233  for (int i = 0; i < ValV.Len(); i++) { NewValV[i] = ValV[i]; }
234  return Entropy(NewValV);
235 }
236 
237 // Entropy of a distribution: ValV[i] contains the number of events i
238 double TSpecFunc::Entropy(const TFltV& ValV) {
239  double Sum=0, Ent=0;
240  for (int i = 0; i < ValV.Len(); i++) {
241  const double& Val = ValV[i];
242  if (Val > 0.0) { Ent -= Val * log(Val); Sum += Val; }
243  }
244  if (Sum > 0.0) {
245  Ent /= Sum;
246  Ent += log(Sum);
247  Ent /= TMath::LogOf2;
248  } else return 1.0;
249  return Ent;
250 }
251 
252 void TSpecFunc::EntropyFracDim(const TIntV& ValV, TFltV& EntropyV) {
253  TFltV NewValV(ValV.Len());
254  for (int i = 0; i < ValV.Len(); i++) {
255  IAssert(ValV[i]==1 || ValV[i] == 0);
256  NewValV[i] = ValV[i];
257  }
258  EntropyFracDim(NewValV, EntropyV);
259 }
260 
261 // Entropy fractal dimension. Input is a vector {0,1}^n.
262 // Where 0 means the event did not occur, and 1 means it occured.
263 // Works exactly as Mengzi Wang's code.
264 void TSpecFunc::EntropyFracDim(const TFltV& ValV, TFltV& EntropyV) {
265  TFltV ValV1, ValV2;
266  int Pow2 = 1;
267  while (2*Pow2 <= ValV.Len()) { Pow2 *= 2; }
268  ValV1.Gen(Pow2);
269  for (int i = 0; i < Pow2; i++) { ValV1[i] = ValV[i];
270  IAssert(ValV[i]==1.0 || ValV[i] == 0.0); }
271  EntropyV.Clr();
272  EntropyV.Add(Entropy(ValV1)); // 2^Pow2 windows
273  while (ValV1.Len() > 2) {
274  ValV2.Gen(ValV1.Len() / 2);
275  for (int i = 0; i < ValV1.Len(); i++) {
276  ValV2[i/2] += ValV1[i];
277  }
278  EntropyV.Add(Entropy(ValV2));
279  ValV1.MoveFrom(ValV2);
280  }
281  EntropyV.Reverse();
282 }
283 
284 // solves for p: B = p * log2(p) + (1-p) * log2(1-p)
285 double TSpecFunc::EntropyBias(const double& B){
286  static TFltFltH BToP;
287  if (BToP.Empty()) {
288  for (double p = 0.5; p < 1.0; p +=0.0001) {
289  double H = p * log(p) + (1.0-p) * log(1.0 - p);
290  H = -H / log(2.0);
291  BToP.AddDat(TMath::Round(H, 3), p);
292  }
293  }
294  if (BToP.IsKey(TMath::Round(B, 3))) { return BToP.GetDat(TMath::Round(B, 3)); }
295  else { return -1.0; }
296 }
297 
298 // MLE of the power-law coefficient
299 double TSpecFunc::GetPowerCoef(const TFltV& XValV, double MinX) {
300  for (int i = 0; MinX <= 0.0 && i < XValV.Len(); i++) {
301  MinX = XValV[i]; }
302  IAssert(MinX > 0.0);
303  double LnSum=0.0;
304  for (int i = 0; i < XValV.Len(); i++) {
305  if (XValV[i].Val < MinX) continue;
306  LnSum += log(XValV[i] / MinX);
307  }
308  return 1.0 + double(XValV.Len()) / LnSum;
309 }
310 
311 double TSpecFunc::GetPowerCoef(const TFltPrV& XValCntV, double MinX) {
312  for (int i = 0; MinX <= 0.0 && i < XValCntV.Len(); i++) {
313  MinX = XValCntV[i].Val1; }
314  IAssert(MinX > 0.0);
315  double NSamples=0.0, LnSum=0.0;
316  for (int i = 0; i < XValCntV.Len(); i++) {
317  if (XValCntV[i].Val1() < MinX) continue;
318  LnSum += XValCntV[i].Val2 * log(XValCntV[i].Val1 / MinX);
319  NSamples += XValCntV[i].Val2;
320  }
321  return 1.0 + NSamples / LnSum;
322 }
323 
325 // Statistical-Moments
326 TMom::TMom(const TFltV& _ValV):
327  //WgtV(_ValV.Len(), 0), ValV(_ValV.Len(), 0),
328  ValWgtV(_ValV.Len(), 0),
329  SumW(), ValSumW(),
330  UsableP(false), UnusableVal(-1),
331  Mn(), Mx(),
332  Mean(), Vari(), SDev(), SErr(),
333  Median(), Quart1(), Quart3(),
334  DecileV(), PercentileV(){
335  for (int ValN=0; ValN<_ValV.Len(); ValN++){Add(_ValV[ValN], 1);}
336  Def();
337 }
338 
339 void TMom::Def(){
340  IAssert(!DefP); DefP=true;
341  UsableP=(SumW>0)&&(ValWgtV.Len()>0);
342  if (UsableP){
343  // Mn, Mx
344  Mn=ValWgtV[0].Val1;
345  Mx=ValWgtV[0].Val1;
346  // Mean, Variance (Mn, Mx), Standard-Error
347  Mean=ValSumW/SumW;
348  Vari=0;
349  if (ValWgtV.Len()>1){
350  for (int ValN=0; ValN<ValWgtV.Len(); ValN++){
351  const double Val=ValWgtV[ValN].Val1;
352  Vari+=ValWgtV[ValN].Val2*TMath::Sqr(Val-Mean);
353  if (Val<Mn){Mn=Val;}
354  if (Val>Mx){Mx=Val;}
355  }
356  Vari=Vari/SumW;
357  // SErr=sqrt(Vari/(ValWgtV.Len()*(ValWgtV.Len()-1))); //J: This seems to be wrong
358  if (Vari > 0.0 && SumW > 0.0) {
359  SErr=sqrt(double(Vari))/sqrt(double(SumW)); //J: This seems to ok: StdDev/sqrt(N)
360  } else { SErr = Mx; } // some big number
361  }
362  // Standard-Deviation
363  SDev=sqrt(double(Vari));
364  // Median
365  ValWgtV.Sort();
366  double CurSumW = 0;
367  for (int ValN=0; ValN<ValWgtV.Len(); ValN++){
368  CurSumW += ValWgtV[ValN].Val2;
369  if (CurSumW > 0.5*SumW) {
370  Median = ValWgtV[ValN].Val1; break; }
371  else if (CurSumW == 0.5*SumW) {
372  Median = 0.5 * (ValWgtV[ValN].Val1+ValWgtV[ValN+1].Val1); break; }
373  }
374  // Quartile-1 and Quartile-3
376  CurSumW = 0;
377  for (int ValN=0; ValN<ValWgtV.Len(); ValN++){
378  CurSumW += ValWgtV[ValN].Val2;
379  if (Quart1==TFlt::Mn) {
380  if (CurSumW > 0.25*SumW) { Quart1 = ValWgtV[ValN].Val1; }
381  //else if (CurSumW == 0.25*SumW) { Quart1 = 0.5 * (ValWgtV[ValN].Val1+ValWgtV[ValN+1].Val1); }
382  }
383  if (Quart3==TFlt::Mn) {
384  if (CurSumW > 0.75*SumW) { Quart3 = ValWgtV[ValN].Val1; }
385  //else if (CurSumW == 0.75*SumW) { Quart3 = 0.5 * (ValWgtV[ValN].Val1+ValWgtV[ValN+1].Val1); }
386  }
387  }
388  // Mode (value with max total weight)
389  THash<TFlt, TFlt> ValWgtH;
390  for (int i = 0; i < ValWgtV.Len(); i++) {
391  ValWgtH.AddDat(ValWgtV[i].Val1) += ValWgtV[i].Val2; }
392  Mode = TFlt::Mn; double MxWgt=TFlt::Mn;
393  for (int v = 0; v < ValWgtH.Len(); v++) {
394  if (ValWgtH[v] > MxWgt) { MxWgt=ValWgtH[v]; Mode=ValWgtH.GetKey(v); }
395  }
396  // Deciles & Percentiles
397  CurSumW = 0;
398  int DecileN = 1, PercentileN = 1;
399  DecileV.Gen(11); PercentileV.Gen(101);
400  DecileV[0]=Mn; DecileV[10]=Mx;
401  PercentileV[0]=Mn; PercentileV[100]=Mx;
402  for (int ValN=0; ValN<ValWgtV.Len(); ValN++){
403  CurSumW += ValWgtV[ValN].Val2;
404  if (CurSumW > SumW*DecileN*0.1) {
405  DecileV[DecileN] = ValWgtV[ValN].Val1; DecileN++; }
406  if (CurSumW > SumW*PercentileN*0.01) {
407  PercentileV[PercentileN] = ValWgtV[ValN].Val1; PercentileN++; }
408  }
409  }
410  ValWgtV.Clr();
411 }
412 
413 
414 
415 double TMom::GetByNm(const TStr& MomNm) const {
416  if (MomNm=="Mean"){return GetMean();}
417  else if (MomNm=="Vari"){return GetVari();}
418  else if (MomNm=="SDev"){return GetSDev();}
419  else if (MomNm=="SErr"){return GetSErr();}
420  else if (MomNm=="Median"){return GetMedian();}
421  else if (MomNm=="Quart1"){return GetQuart1();}
422  else if (MomNm=="Quart3"){return GetQuart3();}
423  else if (MomNm=="Decile0"){return GetDecile(0);}
424  else if (MomNm=="Decile1"){return GetDecile(1);}
425  else if (MomNm=="Decile2"){return GetDecile(2);}
426  else if (MomNm=="Decile3"){return GetDecile(3);}
427  else if (MomNm=="Decile4"){return GetDecile(4);}
428  else if (MomNm=="Decile5"){return GetDecile(5);}
429  else if (MomNm=="Decile6"){return GetDecile(6);}
430  else if (MomNm=="Decile7"){return GetDecile(7);}
431  else if (MomNm=="Decile8"){return GetDecile(8);}
432  else if (MomNm=="Decile9"){return GetDecile(9);}
433  else if (MomNm=="Decile10"){return GetDecile(10);}
434  else {Fail; return 0;}
435 }
436 
437 TStr TMom::GetStrByNm(const TStr& MomNm, char* FmtStr) const {
438  if (IsUsable()){
439  if (FmtStr==NULL){
440  return TFlt::GetStr(GetByNm(MomNm));
441  } else {
442  return TFlt::GetStr(GetByNm(MomNm), FmtStr);
443  }
444  } else {
445  return "X";
446  }
447 }
448 
450  const char& SepCh, const char& DelimCh,
451  const bool& DecileP, const bool& PercentileP, const TStr& FmtStr) const {
452  TChA ChA;
453  if (IsUsable()){
454  ChA+="["; ChA+=SepCh;
455  ChA+="Vals"; ChA+=DelimCh; ChA+=TInt::GetStr(GetVals()); ChA+=SepCh;
456  ChA+="Min"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetMn(), FmtStr.CStr()); ChA+=SepCh;
457  ChA+="Max"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetMx(), FmtStr.CStr()); ChA+=SepCh;
458  ChA+="Mean"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetMean(), FmtStr.CStr()); ChA+=SepCh;
459  //ChA+="Vari"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetVari(), FmtStr.CStr()); ChA+=SepCh;
460  ChA+="SDev"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetSDev(), FmtStr.CStr()); ChA+=SepCh;
461  //ChA+="SErr"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetSErr(), FmtStr.CStr()); ChA+=SepCh;
462  ChA+="Quart1"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetQuart1(), FmtStr.CStr()); ChA+=SepCh;
463  ChA+="Median"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetMedian(), FmtStr.CStr()); ChA+=SepCh;
464  ChA+="Quart3"; ChA+=DelimCh; ChA+=TFlt::GetStr(GetQuart3(), FmtStr.CStr()); ChA+=SepCh;
465  if (DecileP){
466  for (int DecileN=0; DecileN<=10; DecileN++){
467  ChA+="Dec"; ChA+=TInt::GetStr(DecileN);
468  ChA+=DelimCh; ChA+=TFlt::GetStr(GetDecile(DecileN), FmtStr.CStr());
469  ChA+=SepCh;
470  }
471  }
472  if (PercentileP){
473  for (int PercentileN=0; PercentileN<=100; PercentileN++){
474  ChA+="Per"; ChA+=TInt::GetStr(PercentileN);
475  ChA+=DelimCh; ChA+=TFlt::GetStr(GetPercentile(PercentileN), FmtStr.CStr());
476  ChA+=SepCh;
477  }
478  }
479  ChA+="]";
480  } else {
481  ChA="[Unusable]";
482  }
483  return ChA;
484 }
485 
486 TStr TMom::GetNmVStr(const TStr& VarPfx,
487  const char& SepCh, const bool& DecileP, const bool& PercentileP){
488  TChA ChA;
489  ChA+=VarPfx; ChA+="Vals"; ChA+=SepCh;
490  ChA+=VarPfx; ChA+="Min"; ChA+=SepCh;
491  ChA+=VarPfx; ChA+="Max"; ChA+=SepCh;
492  ChA+=VarPfx; ChA+="Mean"; ChA+=SepCh;
493  //ChA+=VarPfx; ChA+="Vari"; ChA+=SepCh;
494  ChA+=VarPfx; ChA+="SDev"; ChA+=SepCh;
495  //ChA+=VarPfx; ChA+="SErr"; ChA+=SepCh;
496  ChA+=VarPfx; ChA+="Quart1"; ChA+=SepCh;
497  ChA+=VarPfx; ChA+="Median"; ChA+=SepCh;
498  ChA+=VarPfx; ChA+="Quart3";
499  if (DecileP){
500  ChA+=SepCh;
501  for (int DecileN=0; DecileN<=10; DecileN++){
502  ChA+=VarPfx; ChA+="Dec"; ChA+=TInt::GetStr(DecileN);
503  if (DecileN<10){ChA+=SepCh;}
504  }
505  }
506  if (PercentileP){
507  ChA+=SepCh;
508  for (int PercentileN=0; PercentileN<=100; PercentileN++){
509  ChA+=VarPfx; ChA+="Per"; ChA+=TInt::GetStr(PercentileN);
510  if (PercentileN<100){ChA+=SepCh;}
511  }
512  }
513  return ChA;
514 }
515 
517  const char& SepCh, const bool& DecileP, const bool& PercentileP) const {
518  TChA ChA;
519  if (IsUsable()){
520  ChA+=TInt::GetStr(GetVals()); ChA+=SepCh;
521  ChA+=TFlt::GetStr(GetMn()); ChA+=SepCh;
522  ChA+=TFlt::GetStr(GetMx()); ChA+=SepCh;
523  ChA+=TFlt::GetStr(GetMean()); ChA+=SepCh;
524  //ChA+=TFlt::GetStr(GetVari()); ChA+=SepCh;
525  ChA+=TFlt::GetStr(GetSDev()); ChA+=SepCh;
526  //ChA+=TFlt::GetStr(GetSErr()); ChA+=SepCh;
527  ChA+=TFlt::GetStr(GetQuart1()); ChA+=SepCh;
528  ChA+=TFlt::GetStr(GetMedian()); ChA+=SepCh;
529  ChA+=TFlt::GetStr(GetQuart3()); ChA+=SepCh;
530  if (DecileP){
531  for (int DecileN=0; DecileN<=10; DecileN++){
532  ChA+=TFlt::GetStr(GetDecile(DecileN)); ChA+=SepCh;
533  }
534  }
535  if (PercentileP){
536  for (int PercentileN=0; PercentileN<=100; PercentileN++){
537  ChA+=TFlt::GetStr(GetPercentile(PercentileN)); ChA+=SepCh;
538  }
539  }
540  } else {
541  int Vals=8;
542  if (DecileP){Vals+=11;}
543  if (PercentileP){Vals+=101;}
544  for (int ValN=0; ValN<Vals; ValN++){
545  ChA="[Unusable]";
546  if (ValN<Vals-1){ChA+=SepCh;}
547  }
548  }
549  return ChA;
550 }
551 
553 // Correlation
554 TCorr::TCorr(const TFltV& ValV1, const TFltV& ValV2):
555  ValVLen(ValV1.Len()), CorrCf(), CorrCfPrb(), FisherZ(){
556  static const double TINY=1.0e-20;
557  IAssert(ValV1.Len()==ValV2.Len());
558 
559  // calculate the means
560  double MeanVal1=0; double MeanVal2=0;
561  {for (int ValN=0; ValN<ValVLen; ValN++){
562  MeanVal1+=ValV1[ValN]; MeanVal2+=ValV2[ValN];}}
563  MeanVal1/=ValVLen; MeanVal2/=ValVLen;
564 
565  // calculate correlation coefficient
566  double yt, xt;
567  double syy=0.0; double sxy=0.0; double sxx=0.0;
568  {for (int ValN=0; ValN<ValVLen; ValN++){
569  xt=ValV1[ValN]-MeanVal1;
570  yt=ValV2[ValN]-MeanVal2;
571  sxx+=xt*xt;
572  syy+=yt*yt;
573  sxy+=xt*yt;
574  }}
575  if (sxx*syy==0){
576  CorrCf=0; //** not in numerical recipes - check why (pojavi se, ko so same nicle)
577  } else {
578  CorrCf=sxy/sqrt(sxx*syy);
579  }
580  // calculate correlation coefficient significance level
581  double df=ValVLen-2;
582  double t=CorrCf*sqrt(df/((1.0-CorrCf+TINY)*(1.0+CorrCf+TINY)));
583  CorrCfPrb=TSpecFunc::BetaI(0.5*df,0.5,df/(df+t*t));
584  // calculate Fisher's Z transformation
585  FisherZ=0.5*log((1.0+(CorrCf)+TINY)/(1.0-(CorrCf)+TINY));
586 }
587 
589 // Statistical Tests
590 void TStatTest::AveVar(const TFltV& ValV, double& Ave, double& Var){
591  Ave=0;
592  for (int ValN=0; ValN<ValV.Len(); ValN++){
593  Ave+=ValV[ValN];}
594  Ave/=ValV.Len();
595  Var=0;
596  double ep=0;
597  for (int ValN=0; ValN<ValV.Len(); ValN++){
598  double s=ValV[ValN]-Ave;
599  ep+=s;
600  Var+=s*s;
601  }
602  Var=(Var-ep*ep/ValV.Len())/(ValV.Len()-1);
603 }
604 
605 double TStatTest::KsProb(const double& Alam) {
606  const double EPS1 = 0.001;
607  const double EPS2 = 1.0e-8;
608  double a2 = -2.0*Alam*Alam, fac = 2.0, sum = 0.0, term, termbf = 0.0;
609  for (int j=1; j <= 100; j++) {
610  term = fac*exp(a2*j*j);
611  sum += term;
612  if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum)
613  return sum;
614  fac = -fac;
615  termbf = fabs(term);
616  }
617  return 1.0;
618 }
619 
621  const TFltV& ObservedBinV, const TFltV& ExpectedBinV,
622  double& ChiSquareVal, double& SignificancePrb){
623  IAssert(ObservedBinV.Len()==ExpectedBinV.Len());
624  int Bins=ObservedBinV.Len();
625  int Constraints=0;
626  int DegreesOfFreedom=Bins-Constraints;
627  ChiSquareVal=0.0;
628  for (int BinN=0; BinN<Bins; BinN++){
629  IAssert(ExpectedBinV[BinN]>0);
630  double BinDiff=ObservedBinV[BinN]-ExpectedBinV[BinN];
631  ChiSquareVal+=BinDiff*BinDiff/ExpectedBinV[BinN];
632  }
633  SignificancePrb=
634  TSpecFunc::GammaQ(0.5*(DegreesOfFreedom), 0.5*(ChiSquareVal));
635 }
636 
638  const TFltV& ObservedBin1V, const TFltV& ObservedBin2V,
639  double& ChiSquareVal, double& SignificancePrb){
640  IAssert(ObservedBin1V.Len()==ObservedBin1V.Len());
641  int Bins=ObservedBin1V.Len();
642  int Constraints=0;
643  int DegreesOfFreedom=Bins-Constraints;
644  ChiSquareVal=0.0;
645  for (int BinN=0; BinN<Bins; BinN++){
646  if ((ObservedBin1V[BinN]==0.0) && (ObservedBin2V[BinN]==0.0)){
647  DegreesOfFreedom--;
648  } else {
649  double BinDiff=ObservedBin1V[BinN]-ObservedBin2V[BinN];
650  ChiSquareVal+=BinDiff*BinDiff/(ObservedBin1V[BinN]+ObservedBin2V[BinN]);
651  }
652  }
653  SignificancePrb=
654  TSpecFunc::GammaQ(0.5*(DegreesOfFreedom),0.5*(ChiSquareVal));
655 }
656 
658  const TFltV& ValV1, const TFltV& ValV2, double& TTestVal, double& TTestPrb){
659  /*double Ave1; double Var1;
660  AveVar(ValV1, Ave1, Var1);
661  double Ave2; double Var2;
662  AveVar(ValV2, Ave2, Var2);*/
663 
664  PMom Val1Mom=TMom::New(ValV1);
665  PMom Val2Mom=TMom::New(ValV2);
666  double ave1=Val1Mom->GetMean();
667  double ave2=Val2Mom->GetMean();
668  double var1=Val1Mom->GetVari();
669  double var2=Val2Mom->GetVari();
670  int n1=ValV1.Len();
671  int n2=ValV2.Len();
672 
673  TTestVal=(ave1-ave2)/sqrt(var1/n1+var2/n2);
674  double df=TMath::Sqr(var1/n1+var2/n2)/(TMath::Sqr(var1/n1)/(n1-1)+TMath::Sqr(var2/n2)/(n2-1));
675  TTestPrb=TSpecFunc::BetaI(0.5*df, 0.5, df/(df+TMath::Sqr(TTestVal)));
676 }
677 
678 void TStatTest::KsTest(const TFltV& ValV1, const TFltV& ValV2, double& DStat, double& PVal) {
679  IAssert(ValV1.IsSorted() && ValV2.IsSorted());
680  int i1=0, i2=0;
681  double CumSum1=0.0, CumSum2=0.0, Cdf1=0.0, Cdf2=0.0;
682  const double N1 = ValV1.Len();
683  const double N2 = ValV2.Len();
684  if (! (N1 > 0.0 && N2 > 0.0)) { DStat = 1.0; PVal = 0.0; return; }
685  DStat=0.0; PVal=0.0;
686  while (i1 < ValV1.Len() && i2 < ValV2.Len()) {
687  const double X1 = ValV1[i1];
688  const double X2 = ValV2[i2];
689  if (X1 <= X2) {
690  CumSum1 += 1;
691  Cdf1 = (CumSum1 / N1);
692  i1++;
693  }
694  if (X2 <= X1) {
695  CumSum2 += 1;
696  Cdf2 = (CumSum2 / N2);
697  i2++;
698  }
699  DStat = TMath::Mx(DStat, fabs(Cdf1 - Cdf2));
700  }
701  const double En = sqrt( N1*N2 / (N1+N2));
702  PVal = TStatTest::KsProb((En+0.12+0.11/En)*DStat);
703 }
704 
705 void TStatTest::KsTest(const TFltPrV& ValCntV1, const TFltPrV& ValCntV2, double& DStat, double& PVal) {
706  IAssert(ValCntV1.IsSorted() && ValCntV2.IsSorted());
707  int i1=0, i2=0;
708  double N1=0.0, N2=0.0, CumSum1=0.0, CumSum2=0.0, Cdf1=0.0, Cdf2=0.0;
709  DStat=0.0; PVal=0.0;
710  for (int i = 0; i < ValCntV1.Len(); i++) N1 += ValCntV1[i].Val2;
711  for (int i = 0; i < ValCntV2.Len(); i++) N2 += ValCntV2[i].Val2;
712  if (! (N1 > 0.0 && N2 > 0.0)) { DStat = 1.0; PVal = 0.0; return; }
713 
714  while (i1 < ValCntV1.Len() && i2 < ValCntV2.Len()) {
715  const double X1 = ValCntV1[i1].Val1;
716  const double X2 = ValCntV2[i2].Val1;
717  if (X1 <= X2) {
718  CumSum1 += ValCntV1[i1].Val2;
719  Cdf1 = (CumSum1 / N1);
720  i1++;
721  }
722  if (X2 <= X1) {
723  CumSum2 += ValCntV2[i2].Val2;
724  Cdf2 = (CumSum2 / N2);
725  i2++;
726  }
727  DStat = TMath::Mx(DStat, fabs(Cdf1 - Cdf2));
728  }
729  const double En = sqrt( N1*N2 / (N1+N2));
730  PVal = TStatTest::KsProb((En+0.12+0.11/En)*DStat);
731 }
732 
734 // Combinations
736  if (ItemV.Len()==0){
737  ItemV.Gen(Order, Order);
738  for (int OrderN=0; OrderN<Order; OrderN++){
739  ItemV[OrderN]=OrderN;}
740  return true;
741  } else {
742  if (ItemV.Last()==Items-1){
743  int OrderN=Order-1;
744  while ((OrderN>=0)&&(ItemV[OrderN]==Items-(Order-OrderN-1)-1)){OrderN--;}
745  if (OrderN<0){
746  return false;
747  } else {
748  ItemV[OrderN]++;
749  for (int SubOrderN=OrderN+1; SubOrderN<Order; SubOrderN++){
750  ItemV[SubOrderN]=ItemV[SubOrderN-1]+1;}
751  CombN++; return true;
752  }
753  } else {
754  ItemV.Last()++; CombN++; return true;
755  }
756  }
757 }
758 
759 int TComb::GetCombs() const {
760  int LCombs=1; int HCombs=1;
761  for (int OrderN=0; OrderN<Order; OrderN++){
762  LCombs*=OrderN+1; HCombs*=Items-OrderN;}
763  int Combs=HCombs/LCombs;
764  return Combs;
765 }
766 
767 void TComb::Wr(){
768  printf("%d:[", GetCombN());
769  for (int OrderN=0; OrderN<Order; OrderN++){
770  if (OrderN>0){printf(" ");}
771  printf("%d", ItemV[OrderN]());
772  }
773  printf("]\n");
774 }
775 
777 // Linear-Regression
778 PLinReg TLinReg::New(const TFltVV& _XVV, const TFltV& _YV, const TFltV& _SigV){
779  PLinReg LinReg=PLinReg(new TLinReg());
780  LinReg->XVV=_XVV;
781  LinReg->YV=_YV;
782  if (_SigV.Empty()){
783  LinReg->SigV.Gen(LinReg->YV.Len());
784  LinReg->SigV.PutAll(1);
785  } else {
786  LinReg->SigV=_SigV;
787  }
788  LinReg->Recs=LinReg->XVV.GetXDim();
789  LinReg->Vars=LinReg->XVV.GetYDim();
790  IAssert(LinReg->Recs>0);
791  IAssert(LinReg->Vars>0);
792  IAssert(LinReg->YV.Len()==LinReg->Recs);
793  IAssert(LinReg->SigV.Len()==LinReg->Recs);
794  LinReg->CovarVV.Gen(LinReg->Vars+1, LinReg->Vars+1);
795  LinReg->CfV.Gen(LinReg->Vars+1);
796  LinReg->NR_lfit();
797  return LinReg;
798 }
799 
801  TFltVV& CovarVV, const int& Vars, const TIntV& ia, const int& mfit){
802  for (int i=mfit+1; i<=Vars; i++){
803  for (int j=1; j<=i; j++){
804  CovarVV.At(i, j)=0; CovarVV.At(j, i)=0.0;}
805  }
806  int k=mfit;
807  for (int j=Vars; j>=1; j--){
808  if (ia[j]!=0){
809  for (int i=1; i<=Vars; i++){Swap(CovarVV.At(i, k), CovarVV.At(i, j));}
810  {for (int i=1; i<=Vars; i++){Swap(CovarVV.At(k, i), CovarVV.At(j, i));}}
811  k--;
812  }
813  }
814 }
815 
816 void TLinReg::NR_gaussj(TFltVV& a, const int& n, TFltVV& b, const int& m){
817  int i, icol=0, irow=0, j, k, l, ll;
818  double big, dum, pivinv;
819 
820  TIntV indxc(n+1);
821  TIntV indxr(n+1);
822  TIntV ipiv(n+1);
823  for (j=1; j<=n; j++){ipiv[j]=0;}
824  for (i=1; i<=n; i++){
825  big=0.0;
826  for (j=1; j<=n; j++){
827  if (ipiv[j]!=1){
828  for (k=1; k<=n; k++){
829  if (ipiv[k]==0){
830  if (fabs(double(a.At(j, k))) >= big){
831  big=fabs(double(a.At(j, k)));
832  irow=j;
833  icol=k;
834  }
835  } else
836  if (ipiv[k]>1){
837  TExcept::Throw("Singular Matrix(1) in Gauss");}
838  }
839  }
840  }
841  ipiv[icol]++;
842  if (irow != icol){
843  for (l=1; l<=n; l++){Swap(a.At(irow, l), a.At(icol, l));}
844  for (l=1; l<=m; l++){Swap(b.At(irow, l), b.At(icol, l));}
845  }
846  indxr[i]=irow;
847  indxc[i]=icol;
848  if (a.At(icol, icol)==0.0){
849  TExcept::Throw("Singular Matrix(1) in Gauss");}
850  pivinv=1.0/a.At(icol, icol);
851  a.At(icol, icol)=1.0;
852  for (l=1; l<=n; l++){a.At(icol, l)=a.At(icol, l)*pivinv;}
853  for (l=1; l<=m; l++){b.At(icol, l)=b.At(icol, l)*pivinv;}
854  for (ll=1; ll<=n; ll++){
855  if (ll != icol){
856  dum=a.At(ll, icol);
857  a.At(ll, icol)=0.0;
858  for (l=1;l<=n;l++){a.At(ll, l)-=a.At(icol, l)*dum;}
859  for (l=1;l<=m;l++){b.At(ll, l)-=b.At(icol, l)*dum;}
860  }
861  }
862  }
863  for (l=n; l>=1; l--){
864  if (indxr[l]!=indxc[l]){
865  for (k=1; k<=n; k++){
866  Swap(a.At(k, indxr[l]), a.At(k, indxc[l]));}
867  }
868  }
869 }
870 
872  int i,j,k,l,m,mfit=0;
873  double ym,wt,sum,sig2i;
874 
875  TIntV ia(Vars+1); for (i=1; i<=Vars; i++){ia[i]=1;}
876  TFltVV beta(Vars+1, 1+1);
877  TFltV afunc(Vars+1);
878  for (j=1;j<=Vars;j++){
879  if (ia[j]!=0){mfit++;}}
880  if (mfit==0){TExcept::Throw("No parameters to be fitted in LFit");}
881  for (j=1; j<=mfit; j++){
882  for (k=1; k<=mfit; k++){CovarVV.At(j, k)=0.0;}
883  beta.At(j, 1)=0.0;
884  }
885  for (i=1; i<=Recs; i++){
886  GetXV(i, afunc); // funcs(XVV[i],afunc,Vars);
887  ym=GetY(i);
888  if (mfit<Vars){
889  for (j=1;j<=Vars;j++){
890  if (ia[j]==0){ym-=CfV[j]*afunc[j];}}
891  }
892  sig2i=1.0/TMath::Sqr(GetSig(i));
893  for (j=0, l=1; l<=Vars; l++){
894  if (ia[l]!=0){
895  wt=afunc[l]*sig2i;
896  for (j++, k=0, m=1; m<=l; m++){
897  if (ia[m]!=0){CovarVV.At(j, ++k)+=wt*afunc[m];}
898  }
899  beta.At(j, 1)+=ym*wt;
900  }
901  }
902  }
903  for (j=2; j<=mfit; j++){
904  for (k=1; k<j; k++){CovarVV.At(k, j)=CovarVV.At(j, k);}
905  }
906  NR_gaussj(CovarVV, mfit, beta, 1);
907  for (j=0, l=1; l<=Vars; l++){
908  if (ia[l]!=0){CfV[l]=beta.At(++j, 1);}
909  }
910  ChiSq=0.0;
911  for (i=1; i<=Recs; i++){
912  GetXV(i, afunc); // funcs(XVV[i],afunc,Vars);
913  for (sum=0.0, j=1; j<=Vars; j++){sum+=CfV[j]*afunc[j];}
914  ChiSq+=TMath::Sqr((GetY(i)-sum)/GetSig(i));
915  }
916  NR_covsrt(CovarVV, Vars, ia, mfit);
917 }
918 
919 void TLinReg::Wr() const {
920  printf("\n%11s %21s\n","parameter","uncertainty");
921  for (int i=0; i<Vars;i++){
922  printf(" a[%1d] = %8.6f %12.6f\n",
923  i+1, GetCf(i), GetCfUncer(i));
924  }
925  printf("chi-squared = %12f\n", GetChiSq());
926 
927  printf("full covariance matrix\n");
928  {for (int i=0;i<Vars;i++){
929  for (int j=0;j<Vars;j++){
930  printf("%12f", GetCovar(i, j));}
931  printf("\n");
932  }}
933 }
934 
936 // Singular-Value-Decomposition
937 PSvd TSvd::New(const TFltVV& _XVV, const TFltV& _YV, const TFltV& _SigV){
938  PSvd Svd=PSvd(new TSvd());
939  Svd->XVV=_XVV;
940  Svd->YV=_YV;
941  if (_SigV.Empty()){
942  Svd->SigV.Gen(Svd->YV.Len());
943  Svd->SigV.PutAll(1);
944  } else {
945  Svd->SigV=_SigV;
946  }
947  Svd->Recs=Svd->XVV.GetXDim();
948  Svd->Vars=Svd->XVV.GetYDim();
949  IAssert(Svd->Recs>0);
950  IAssert(Svd->Vars>0);
951  IAssert(Svd->YV.Len()==Svd->Recs);
952  IAssert(Svd->SigV.Len()==Svd->Recs);
953  Svd->CovarVV.Gen(Svd->Vars+1, Svd->Vars+1);
954  Svd->CfV.Gen(Svd->Vars+1);
955  Svd->NR_svdfit();
956  return Svd;
957 }
958 
959 double TSvd::NR_pythag(double a, double b){
960  double absa,absb;
961  absa=fabs(a);
962  absb=fabs(b);
963  if (absa > absb){
964  return absa*sqrt(1.0+TMath::Sqr(absb/absa));
965  } else {
966  return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+TMath::Sqr(absa/absb)));
967  }
968 }
969 
970 void TSvd::NR_svdcmp(TFltVV& a, int m, int n, TFltV& w, TFltVV& v){
971  int flag,i,its,j,jj,k,l=0,nm;
972  double anorm,c,f,g,h,s,scale,x,y,z;
973 
974  TFltV rv1(n+1);
975  g=scale=anorm=0.0;
976  for (i=1;i<=n;i++) {
977  l=i+1;
978  rv1[i]=scale*g;
979  g=s=scale=0.0;
980  if (i <= m) {
981  for (k=i;k<=m;k++) scale += fabs(double(a.At(k,i)));
982  if (scale) {
983  for (k=i;k<=m;k++) {
984  a.At(k,i) /= scale;
985  s += a.At(k,i)*a.At(k,i);
986  }
987  f=a.At(i,i);
988  g = -NR_SIGN(sqrt(s),f);
989  h=f*g-s;
990  a.At(i,i)=f-g;
991  for (j=l;j<=n;j++) {
992  for (s=0.0,k=i;k<=m;k++) s += a.At(k,i)*a(k,j);
993  f=s/h;
994  for (k=i;k<=m;k++) a.At(k,j) += f*a.At(k,i);
995  }
996  for (k=i;k<=m;k++) a.At(k,i) *= scale;
997  }
998  }
999  w[i]=scale *g;
1000  g=s=scale=0.0;
1001  if (i <= m && i != n) {
1002  for (k=l;k<=n;k++) scale += fabs(double(a.At(i,k)));
1003  if (scale) {
1004  for (k=l;k<=n;k++) {
1005  a.At(i,k) /= scale;
1006  s += a.At(i,k)*a.At(i,k);
1007  }
1008  f=a.At(i,l);
1009  g = -NR_SIGN(sqrt(s),f);
1010  h=f*g-s;
1011  a.At(i,l)=f-g;
1012  for (k=l;k<=n;k++) rv1[k]=a.At(i,k)/h;
1013  for (j=l;j<=m;j++) {
1014  for (s=0.0,k=l;k<=n;k++) s += a.At(j,k)*a.At(i,k);
1015  for (k=l;k<=n;k++) a.At(j,k) += s*rv1[k];
1016  }
1017  for (k=l;k<=n;k++) a.At(i,k) *= scale;
1018  }
1019  }
1020  anorm=NR_FMAX(anorm,(fabs(double(w[i]))+fabs(double(rv1[i]))));
1021  }
1022  for (i=n;i>=1;i--) {
1023  if (i < n) {
1024  if (g) {
1025  for (j=l;j<=n;j++)
1026  v.At(j,i)=(a.At(i,j)/a.At(i,l))/g;
1027  for (j=l;j<=n;j++) {
1028  for (s=0.0,k=l;k<=n;k++) s += a.At(i,k)*v.At(k,j);
1029  for (k=l;k<=n;k++) v.At(k,j) += s*v.At(k,i);
1030  }
1031  }
1032  for (j=l;j<=n;j++) v.At(i,j)=v.At(j,i)=0.0;
1033  }
1034  v.At(i,i)=1.0;
1035  g=rv1[i];
1036  l=i;
1037  }
1038  for (i=NR_IMIN(m,n);i>=1;i--) {
1039  l=i+1;
1040  g=w[i];
1041  for (j=l;j<=n;j++) a.At(i,j)=0.0;
1042  if (g) {
1043  g=1.0/g;
1044  for (j=l;j<=n;j++) {
1045  for (s=0.0,k=l;k<=m;k++) s += a.At(k,i)*a.At(k,j);
1046  f=(s/a.At(i,i))*g;
1047  for (k=i;k<=m;k++) a.At(k,j) += f*a.At(k,i);
1048  }
1049  for (j=i;j<=m;j++) a.At(j,i) *= g;
1050  } else for (j=i;j<=m;j++) a.At(j,i)=0.0;
1051  a.At(i,i)++;
1052  }
1053  for (k=n;k>=1;k--) {
1054  for (its=1;its<=30;its++) {
1055  flag=1;
1056  for (l=k;l>=1;l--) {
1057  nm=l-1;
1058  if ((double)(fabs(double(rv1[l])+anorm)) == anorm) {
1059  flag=0;
1060  break;
1061  }
1062  if ((double)(fabs(double(w[nm]))+anorm) == anorm) break;
1063  }
1064  if (flag) {
1065  c=0.0;
1066  s=1.0;
1067  for (i=l;i<=k;i++) {
1068  f=s*rv1[i];
1069  rv1[i]=c*rv1[i];
1070  if ((double)(fabs(f)+anorm) == anorm) break;
1071  g=w[i];
1072  h=NR_pythag(f,g);
1073  w[i]=h;
1074  h=1.0/h;
1075  c=g*h;
1076  s = -f*h;
1077  for (j=1;j<=m;j++) {
1078  y=a.At(j,nm);
1079  z=a.At(j,i);
1080  a.At(j,nm)=y*c+z*s;
1081  a.At(j,i)=z*c-y*s;
1082  }
1083  }
1084  }
1085  z=w[k];
1086  if (l == k) {
1087  if (z < 0.0) {
1088  w[k] = -z;
1089  for (j=1;j<=n;j++) v.At(j,k) = -v.At(j,k);
1090  }
1091  break;
1092  }
1093  if (its==30){
1094  TExcept::Throw("no convergence in 30 svdcmp iterations");}
1095  x=w[l];
1096  nm=k-1;
1097  y=w[nm];
1098  g=rv1[nm];
1099  h=rv1[k];
1100  f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
1101  g=NR_pythag(f,1.0);
1102  f=((x-z)*(x+z)+h*((y/(f+NR_SIGN(g,f)))-h))/x;
1103  c=s=1.0;
1104  for (j=l;j<=nm;j++) {
1105  i=j+1;
1106  g=rv1[i];
1107  y=w[i];
1108  h=s*g;
1109  g=c*g;
1110  z=NR_pythag(f,h);
1111  rv1[j]=z;
1112  c=f/z;
1113  s=h/z;
1114  f=x*c+g*s;
1115  g = g*c-x*s;
1116  h=y*s;
1117  y *= c;
1118  for (jj=1;jj<=n;jj++) {
1119  x=v.At(jj,j);
1120  z=v.At(jj,i);
1121  v.At(jj,j)=x*c+z*s;
1122  v.At(jj,i)=z*c-x*s;
1123  }
1124  z=NR_pythag(f,h);
1125  w[j]=z;
1126  if (z) {
1127  z=1.0/z;
1128  c=f*z;
1129  s=h*z;
1130  }
1131  f=c*g+s*y;
1132  x=c*y-s*g;
1133  for (jj=1;jj<=m;jj++) {
1134  y=a.At(jj,j);
1135  z=a.At(jj,i);
1136  a.At(jj,j)=y*c+z*s;
1137  a.At(jj,i)=z*c-y*s;
1138  }
1139  }
1140  rv1[l]=0.0;
1141  rv1[k]=f;
1142  w[k]=x;
1143  }
1144  }
1145 }
1146 
1148  TFltVV& u, TFltV& w, TFltVV& v, int m, int n, TFltV& b, TFltV& x){
1149  int jj,j,i;
1150  double s;
1151 
1152  TFltV tmp(n+1);
1153  for (j=1;j<=n;j++) {
1154  s=0.0;
1155  if (w[j]) {
1156  for (i=1;i<=m;i++) s += u.At(i,j)*b[i];
1157  s /= w[j];
1158  }
1159  tmp[j]=s;
1160  }
1161  for (j=1;j<=n;j++) {
1162  s=0.0;
1163  for (jj=1;jj<=n;jj++) s += v.At(j,jj)*tmp[jj];
1164  x[j]=s;
1165  }
1166 }
1167 
1168 void TSvd::NR_svdvar(TFltVV& v, int ma, TFltV& w, TFltVV& cvm){
1169  int k,j,i;
1170  double sum;
1171 
1172  TFltV wti(ma+1);
1173  for (i=1;i<=ma;i++) {
1174  wti[i]=0.0;
1175  if (w[i]) wti[i]=1.0/(w[i]*w[i]);
1176  }
1177  for (i=1;i<=ma;i++) {
1178  for (j=1;j<=i;j++) {
1179  for (sum=0.0,k=1;k<=ma;k++) sum += v.At(i,k)*v.At(j,k)*wti[k];
1180  cvm.At(j,i)=cvm.At(i,j)=sum;
1181  }
1182  }
1183 }
1184 
1186  int j,i;
1187  double wmax,tmp,thresh,sum;
1188  double TOL=1.0e-5;
1189 
1190  TFltVV u(Recs+1, Vars+1);
1191  TFltVV v(Vars+1, Vars+1);
1192  TFltV w(Vars+1);
1193  TFltV b(Recs+1);
1194  TFltV afunc(Vars+1);
1195  for (i=1;i<=Recs;i++) {
1196  GetXV(i, afunc); // (*funcs)(x[i],afunc,Vars);
1197  tmp=1.0/GetSig(i);
1198  for (j=1;j<=Vars;j++){u.At(i,j)=afunc[j]*tmp;}
1199  b[i]=GetY(i)*tmp;
1200  }
1201  NR_svdcmp(u,Recs,Vars,w,v);
1202  wmax=0.0;
1203  for (j=1;j<=Vars;j++){
1204  if (w[j] > wmax){wmax=w[j];}}
1205  thresh=TOL*wmax;
1206  for (j=1;j<=Vars;j++){
1207  if (double(w[j])<thresh){w[j]=0.0;}}
1208  NR_svbksb(u,w,v,Recs,Vars,b,CfV);
1209  ChiSq=0.0;
1210  for (i=1;i<=Recs;i++) {
1211  GetXV(i, afunc); // (*funcs)(x[i],afunc,Vars);
1212  for (sum=0.0,j=1;j<=Vars;j++){sum += CfV[j]*afunc[j];}
1213  ChiSq += (tmp=(GetY(i)-sum)/GetSig(i),tmp*tmp);
1214  }
1215  // covariance matrix calculation
1216  CovarVV.Gen(Vars+1, Vars+1);
1217  NR_svdvar(v, Vars, w, CovarVV);
1218 }
1219 
1220 void TSvd::GetCfV(TFltV& _CfV){
1221  _CfV=CfV; _CfV.Del(0);
1222 }
1223 
1224 void TSvd::GetCfUncerV(TFltV& CfUncerV){
1225  CfUncerV.Gen(Vars);
1226  for (int VarN=0; VarN<Vars; VarN++){
1227  CfUncerV[VarN]=GetCfUncer(VarN);
1228  }
1229 }
1230 
1231 // all vectors (matrices) start with 0
1232 void TSvd::Svd(const TFltVV& InMtx, TFltVV& LSingV, TFltV& SingValV, TFltVV& RSingV) {
1233  //LSingV = InMtx;
1234  LSingV.Gen(InMtx.GetYDim()+1, InMtx.GetYDim()+1);
1235  // create 1 based adjacency matrix
1236  for (int x = 0; x < InMtx.GetXDim(); x++) {
1237  for (int y = 0; y < InMtx.GetYDim(); y++) {
1238  LSingV.At(x+1, y+1) = InMtx.At(x, y);
1239  }
1240  }
1241  RSingV.Gen(InMtx.GetYDim()+1, InMtx.GetYDim()+1);
1242  SingValV.Gen(InMtx.GetYDim()+1);
1243  TSvd::NR_svdcmp(LSingV, InMtx.GetXDim(), InMtx.GetYDim(), SingValV, RSingV);
1244  // 0-th singular value/vector is full of zeros, delete it
1245  SingValV.Del(0);
1246  LSingV.DelX(0); LSingV.DelY(0);
1247  RSingV.DelX(0); RSingV.DelY(0);
1248 }
1249 
1250 // InMtx1 is 1-based (0-th row/column are full of zeros)
1251 // returned singular vectors/values are 0 based
1252 void TSvd::Svd1Based(const TFltVV& InMtx1, TFltVV& LSingV, TFltV& SingValV, TFltVV& RSingV) {
1253  LSingV = InMtx1;
1254  SingValV.Gen(InMtx1.GetYDim());
1255  RSingV.Gen(InMtx1.GetYDim(), InMtx1.GetYDim());
1256  TSvd::NR_svdcmp(LSingV, InMtx1.GetXDim()-1, InMtx1.GetYDim()-1, SingValV, RSingV);
1257  // 0-th singular value/vector is full of zeros, delete it
1258  SingValV.Del(0);
1259  LSingV.DelX(0); LSingV.DelY(0);
1260  RSingV.DelX(0); RSingV.DelY(0);
1261 }
1262 
1263 void TSvd::Wr() const {
1264  printf("\n%11s %21s\n","parameter","uncertainty");
1265  for (int i=0; i<Vars;i++){
1266  printf(" a[%1d] = %8.6f %12.6f\n",
1267  i+1, GetCf(i), GetCfUncer(i));
1268  }
1269  printf("chi-squared = %12f\n", GetChiSq());
1270 
1271  printf("full covariance matrix\n");
1272  {for (int i=0;i<Vars;i++){
1273  for (int j=0;j<Vars;j++){
1274  printf("%12f", GetCovar(i, j));}
1275  printf("\n");
1276  }}
1277 }
1278 
1280 // Histogram
1281 void THist::Add(const double& Val, const bool& OnlyInP) {
1282  // get bucket number
1283  const int BucketN = int(floor((Val - MnVal) / BucketSize));
1284  if (OnlyInP) {
1285  // value should be inside
1286  EAssert(MnVal <= Val && Val <= MxVal);
1287  BucketV[BucketN]++;
1288  } else {
1289  // value does not need to be inside
1290  if (BucketN < 0) {
1291  BucketV[0]++;
1292  } else if (BucketN < BucketV.Len()) {
1293  BucketV[BucketN]++;
1294  } else {
1295  BucketV.Last()++;
1296  }
1297  }
1298  // for computing percentage
1299  Vals++;
1300 }
1301 
1302 void THist::SaveStat(const TStr& ValNm, TSOut& FOut) const {
1303  FOut.PutStrLn("#" + ValNm + ": " + Vals.GetStr());
1304  const int Buckets = BucketV.Len() - 1;
1305  for (int BucketN = 0; BucketN < Buckets; BucketN++) {
1306  FOut.PutStrLn(TStr::Fmt("%d-%d\t%d", BucketSize*BucketN,
1307  BucketSize*(BucketN+1), BucketV[BucketN]()));
1308  }
1309  if (BucketV.Last() > 0) {
1310  FOut.PutStrLn(TStr::Fmt("%d-\t%d", BucketSize*Buckets, BucketV.Last()()));
1311  }
1312 
1313 }
1314 
#define IAssert(Cond)
Definition: bd.h:262
TFlt SDev
Definition: xmath.h:138
TBool UsableP
Definition: xmath.h:135
static double LnGamma(const double &xx)
Definition: xmath.cpp:80
double GetCf(const int &VarN) const
Definition: xmath.h:377
TStr GetStr() const
Definition: dt.h:1197
TCorr()
Definition: xmath.h:276
static double NR_pythag(double a, double b)
Definition: xmath.cpp:959
void Wr() const
Definition: xmath.cpp:919
double FisherZ
Definition: xmath.h:274
double GetMedian() const
Definition: xmath.h:244
static TStr GetNmVStr(const TStr &VarPfx, const char &SepCh='\t', const bool &DecileP=true, const bool &PercentileP=true)
Definition: xmath.cpp:486
TFlt MxVal
Definition: xmath.h:459
double GetSig(const int RecN) const
Definition: xmath.h:359
void Del(const TSizeTy &ValN)
Removes the element at position ValN.
Definition: ds.h:1189
double GetCovar(const int &VarN1, const int &VarN2) const
Definition: xmath.h:440
int Vars
Definition: xmath.h:350
static const T & Mx(const T &LVal, const T &RVal)
Definition: xmath.h:32
TInt Vals
Definition: xmath.h:462
static double KsProb(const double &Alam)
Definition: xmath.cpp:605
void NR_gaussj(TFltVV &a, const int &n, TFltVV &b, const int &m)
Definition: xmath.cpp:816
double CorrCf
Definition: xmath.h:272
static void PowerFit(const TVec< TFltPr > &XY, double &A, double &B, double &SigA, double &SigB, double &Chi2, double &R2)
Definition: xmath.cpp:193
int GetVals() const
Definition: xmath.h:220
#define Fail
Definition: bd.h:238
static PLinReg New(const TFltVV &XVV, const TFltV &YV, const TFltV &SigV=TFltV())
Definition: xmath.cpp:778
void Wr()
Definition: xmath.cpp:767
int Recs
Definition: xmath.h:402
static void ExpFit(const TVec< TFltPr > &XY, double &A, double &B, double &SigA, double &SigB, double &Chi2, double &R2)
Definition: xmath.cpp:218
static int NR_IMIN(int iminarg1, int iminarg2)
Definition: xmath.h:415
TSizeTy Len() const
Returns the number of elements in the vector.
Definition: ds.h:575
void NR_svdfit()
Definition: xmath.cpp:1185
TMom()
Definition: xmath.h:144
TFlt Mx
Definition: xmath.h:137
bool Empty() const
Definition: hash.h:227
static void LogFit(const TVec< TFltPr > &XY, double &A, double &B, double &SigA, double &SigB, double &Chi2, double &R2)
Definition: xmath.cpp:208
double GetCovar(const int &VarN1, const int &VarN2) const
Definition: xmath.h:380
int Order
Definition: xmath.h:319
TFlt Mode
Definition: xmath.h:140
void DelX(const TSizeTy &X)
Definition: ds.h:2357
static double NR_FMAX(double maxarg1, double maxarg2)
Definition: xmath.h:413
void GetXV(const int RecN, TFltV &VarV) const
Definition: xmath.h:406
const TDat & GetDat(const TKey &Key) const
Definition: hash.h:262
static double Sqr(const double &x)
Definition: xmath.h:12
void Add(const double &Val, const bool &OnlyInP)
Definition: xmath.cpp:1281
int ValVLen
Definition: xmath.h:271
double GetSDev() const
Definition: xmath.h:242
void SaveStat(const TStr &ValNm, TSOut &FOut) const
Definition: xmath.cpp:1302
double GetChiSq() const
Definition: xmath.h:443
TPt< TSvd > PSvd
Definition: xmath.h:397
bool IsUsable() const
Definition: xmath.h:225
bool Empty() const
Tests whether the vector is empty.
Definition: ds.h:570
static double LogOf2
Definition: xmath.h:9
static double GammaQ(const double &a, const double &x)
Definition: xmath.cpp:68
static void NR_svdcmp(TFltVV &a, int m, int n, TFltV &w, TFltVV &v)
Definition: xmath.cpp:970
int Items
Definition: xmath.h:318
TFlt Mn
Definition: xmath.h:137
double GetSErr() const
Definition: xmath.h:243
static void Svd(const TFltVV &InMtx, TFltVV &LSingV, TFltV &SingValV, TFltVV &RSingV)
Definition: xmath.cpp:1232
void Clr(const bool &DoDel=true, const TSizeTy &NoDelLim=-1)
Clears the contents of the vector.
Definition: ds.h:1022
int Recs
Definition: xmath.h:350
double GetMx() const
Definition: xmath.h:238
TStr GetValVStr(const char &SepCh='\t', const bool &DecileP=true, const bool &PercentileP=true) const
Definition: xmath.cpp:516
void Add(const TFlt &Val, const TFlt &Wgt=1)
Definition: xmath.h:217
double GetCfUncer(const int &VarN) const
Definition: xmath.h:438
TSvd()
Definition: xmath.h:424
static void EntropyFracDim(const TIntV &ValV, TFltV &EntropyV)
Definition: xmath.cpp:252
TFltV CfV
Definition: xmath.h:404
void Sort(const bool &Asc=true)
Sorts the elements of the vector.
Definition: ds.h:1318
void GetCfV(TFltV &_CfV)
Definition: xmath.cpp:1220
void NR_covsrt(TFltVV &CovarVV, const int &Vars, const TIntV &ia, const int &mfit)
Definition: xmath.cpp:800
static void Throw(const TStr &MsgStr)
Definition: ut.h:187
static PMom New()
Definition: xmath.h:160
TSizeTy GetYDim() const
Definition: ds.h:2250
static void ChiSquareTwo(const TFltV &ObservedBin1V, const TFltV &ObservedBin2V, double &ChiSquareVal, double &SignificancePrb)
Definition: xmath.cpp:637
static double Round(const double &Val)
Definition: xmath.h:16
double ChiSq
Definition: xmath.h:353
TFlt Vari
Definition: xmath.h:138
TLinReg()
Definition: xmath.h:364
TFlt SumW
Definition: xmath.h:133
static void Svd1Based(const TFltVV &InMtx1, TFltVV &LSingV, TFltV &SingValV, TFltVV &RSingV)
Definition: xmath.cpp:1252
double ChiSq
Definition: xmath.h:405
TFltV CfV
Definition: xmath.h:352
TFlt ValSumW
Definition: xmath.h:133
const TVal & Last() const
Returns a reference to the last element of the vector.
Definition: ds.h:579
static void ChiSquareOne(const TFltV &ObservedBinV, const TFltV &ExpectedBinV, double &ChiSquareVal, double &SignificancePrb)
Definition: xmath.cpp:620
static void AveVar(const TFltV &ValV, double &Ave, double &Var)
Definition: xmath.cpp:590
void GetCfUncerV(TFltV &CfUncerV)
Definition: xmath.cpp:1224
int GetCombs() const
Definition: xmath.cpp:759
TFltVV CovarVV
Definition: xmath.h:403
TPair< TFlt, TFlt > TFltPr
Definition: ds.h:99
Definition: fl.h:128
static double Pi
Definition: xmath.h:8
TPt< TLinReg > PLinReg
Definition: xmath.h:345
int GetCombN() const
Definition: xmath.h:338
double GetPercentile(const int &PercentileN) const
Definition: xmath.h:250
void NR_svdvar(TFltVV &v, int ma, TFltV &w, TFltVV &cvm)
Definition: xmath.cpp:1168
TFlt Mean
Definition: xmath.h:138
double GetChiSq() const
Definition: xmath.h:383
double GetQuart3() const
Definition: xmath.h:246
TFlt Quart3
Definition: xmath.h:139
static double EntropyBias(const double &B)
Definition: xmath.cpp:285
static void KsTest(const TFltV &ValV1, const TFltV &ValV2, double &DStat, double &PVal)
Definition: xmath.cpp:678
static void GammaQContFrac(double &gammcf, const double &a, const double &x, double &gln)
Definition: xmath.cpp:38
TStr GetStr(const char &SepCh=' ', const char &DelimCh=':', const bool &DecileP=true, const bool &PercentileP=true, const TStr &FmtStr="%g") const
Definition: xmath.cpp:449
static double GetPowerCoef(const TFltV &XValV, double MinX=-1.0)
Definition: xmath.cpp:299
Definition: dt.h:201
#define EAssert(Cond)
Definition: bd.h:280
static double NR_SIGN(double a, double b)
Definition: xmath.h:412
TIntV BucketV
Definition: xmath.h:460
TInt Vals
Definition: xmath.h:134
bool GetNext()
Definition: xmath.cpp:735
static double BetaI(const double &a, const double &b, const double &x)
Definition: xmath.cpp:137
void DelY(const TSizeTy &Y)
Definition: ds.h:2369
double GetMean() const
Definition: xmath.h:240
Definition: dt.h:412
static TStr Fmt(const char *FmtStr,...)
Definition: dt.cpp:1599
bool IsSorted(const bool &Asc=true) const
Checks whether the vector is sorted in ascending (if Asc=true) or descending (if Asc=false) order...
Definition: ds.h:1323
static void TTest(const TFltV &ValV1, const TFltV &ValV2, double &TTestVal, double &TTestPrb)
Definition: xmath.cpp:657
TSizeTy GetXDim() const
Definition: ds.h:2249
double CorrCfPrb
Definition: xmath.h:273
Definition: hash.h:97
static PSvd New(const TFltVV &XVV, const TFltV &YV, const TFltV &SigV=TFltV())
Definition: xmath.cpp:937
static double LnComb(const int &n, const int &k)
Definition: xmath.cpp:95
double GetQuart1() const
Definition: xmath.h:245
void GetXV(const int RecN, TFltV &VarV) const
Definition: xmath.h:354
TStr GetStrByNm(const TStr &MomNm, char *FmtStr=NULL) const
Definition: xmath.cpp:437
static void GammaPSeries(double &gamser, const double &a, const double &x, double &gln)
Definition: xmath.cpp:9
double GetY(const int RecN) const
Definition: xmath.h:358
TFltPrV ValWgtV
Definition: xmath.h:132
double GetVari() const
Definition: xmath.h:241
void NR_svbksb(TFltVV &u, TFltV &w, TFltVV &v, int m, int n, TFltV &b, TFltV &x)
Definition: xmath.cpp:1147
Definition: bd.h:196
static double Entropy(const TIntV &ValV)
Definition: xmath.cpp:231
void Reverse()
Reverses the order of the elements in the vector.
Definition: ds.h:1350
double GetDecile(const int &DecileN) const
Definition: xmath.h:248
double GetCf(const int &VarN) const
Definition: xmath.h:437
static void LinearFit(const TVec< TFltPr > &XY, double &A, double &B, double &SigA, double &SigB, double &Chi2, double &R2)
Definition: xmath.cpp:150
TStr GetStr() const
Definition: dt.h:1459
TFlt SErr
Definition: xmath.h:138
int CombN
Definition: xmath.h:320
void Gen(const TSizeTy &_Vals)
Constructs a vector (an array) of _Vals elements.
Definition: ds.h:523
TFlt BucketSize
Definition: xmath.h:461
double GetByNm(const TStr &MomNm) const
Definition: xmath.cpp:415
void MoveFrom(TVec< TVal, TSizeTy > &Vec)
Takes over the data and the capacity from Vec.
Definition: ds.h:1073
int Vars
Definition: xmath.h:402
TFlt Median
Definition: xmath.h:139
char * CStr()
Definition: dt.h:476
TFltVV CovarVV
Definition: xmath.h:351
bool IsKey(const TKey &Key) const
Definition: hash.h:258
int PutStrLn(const TStr &Str, const bool &ForceInLn=false)
Definition: fl.h:161
void NR_lfit()
Definition: xmath.cpp:871
TBool DefP
Definition: xmath.h:131
TSizeTy Add()
Adds a new element at the end of the vector, after its current last element.
Definition: ds.h:602
double GetY(const int RecN) const
Definition: xmath.h:410
void Wr() const
Definition: xmath.cpp:1263
double GetSig(const int RecN) const
Definition: xmath.h:411
TFltV PercentileV
Definition: xmath.h:142
int Len() const
Definition: hash.h:228
TDat & AddDat(const TKey &Key)
Definition: hash.h:238
static double BetaCf(const double &a, const double &b, const double &x)
Definition: xmath.cpp:99
TFltV DecileV
Definition: xmath.h:141
double GetMn() const
Definition: xmath.h:237
static double E
Definition: xmath.h:7
void Def()
Definition: xmath.cpp:339
void Gen(const TSizeTy &_XDim, const TSizeTy &_YDim)
Definition: ds.h:2246
const TKey & GetKey(const int &KeyId) const
Definition: hash.h:252
TFlt MnVal
Definition: xmath.h:458
TIntV ItemV
Definition: xmath.h:321
void Swap(TRec &Rec1, TRec &Rec2)
Definition: bd.h:568
static const double Mn
Definition: dt.h:1387
TFlt Quart1
Definition: xmath.h:139
const TVal & At(const TSizeTy &X, const TSizeTy &Y) const
Definition: ds.h:2255
double GetCfUncer(const int &VarN) const
Definition: xmath.h:378