SNAP Library 6.0, User Reference  2020-12-09 16:24:20
SNAP, a general purpose, high performance system for analysis and manipulation of large networks
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:1200
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:2358
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:2251
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:2370
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:2250
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:1462
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:479
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:2247
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:1390
TFlt Quart1
Definition: xmath.h:139
const TVal & At(const TSizeTy &X, const TSizeTy &Y) const
Definition: ds.h:2256
double GetCfUncer(const int &VarN) const
Definition: xmath.h:378