SHOGUN  v3.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Statistics.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011-2012 Heiko Strathmann
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  *
10  * ALGLIB Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier under GPL2+
11  * http://www.alglib.net/
12  * See header file for which functions are taken from ALGLIB (with adjustments
13  * for shogun)
14  */
15 
16 #include <shogun/mathematics/Statistics.h>
17 #include <shogun/mathematics/Math.h>
18 #include <shogun/lib/SGMatrix.h>
19 #include <shogun/lib/SGVector.h>
20 #include <shogun/lib/SGSparseMatrix.h>
21 #include <shogun/lib/SGSparseVector.h>
22 
23 #ifdef HAVE_LAPACK
24 #include <shogun/mathematics/lapack.h>
25 #endif //HAVE_LAPACK
26 
27 #ifdef HAVE_EIGEN3
28 #include <shogun/mathematics/eigen3.h>
29 using namespace Eigen;
30 #endif //HAVE_EIGEN3
31 
32 using namespace shogun;
33 
34 float64_t CStatistics::mean(SGVector<float64_t> values)
35 {
36  ASSERT(values.vlen>0)
37  ASSERT(values.vector)
38 
39  float64_t sum=0;
40  for (index_t i=0; i<values.vlen; ++i)
41  sum+=values.vector[i];
42 
43  return sum/values.vlen;
44 }
45 
46 float64_t CStatistics::median(SGVector<float64_t> values, bool modify,
47  bool in_place)
48 {
49  float64_t result;
50  if (modify)
51  {
52  /* use QuickSelect method
53  * This Quickselect routine is based on the algorithm described in
54  * "Numerical recipes in C", Second Edition,
55  * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
56  * This code by Nicolas Devillard - 1998. Public domain.
57  * Adapted to SHOGUN by Heiko Strathmann
58  */
59  int32_t low;
60  int32_t high;
61  int32_t median;
62  int32_t middle;
63  int32_t l;
64  int32_t h;
65 
66  low=0;
67  high=values.vlen-1;
68  median=(low+high)/2;
69 
70  while (true)
71  {
72  if (high<=low)
73  {
74  result=values[median];
75  break;
76  }
77 
78  if (high==low+1)
79  {
80  if (values[low]>values[high])
81  CMath::CMath::swap(values[low], values[high]);
82  result=values[median];
83  break;
84  }
85 
86  middle=(low+high)/2;
87  if (values[middle]>values[high])
88  CMath::swap(values[middle], values[high]);
89  if (values[low]>values[high])
90  CMath::swap(values[low], values[high]);
91  if (values[middle]>values[low])
92  CMath::swap(values[middle], values[low]);
93 
94  CMath::swap(values[middle], values[low+1]);
95 
96  l=low+1;
97  h=high;
98  for (;;)
99  {
100  do
101  l++;
102  while (values[low]>values[l]);
103  do
104  h--;
105  while (values[h]>values[low]);
106  if (h<l)
107  break;
108  CMath::swap(values[l], values[h]);
109  }
110 
111  CMath::swap(values[low], values[h]);
112  if (h<=median)
113  low=l;
114  if (h>=median)
115  high=h-1;
116  }
117 
118  }
119  else
120  {
121  if (in_place)
122  {
123  /* use Torben method
124  * The following code is public domain.
125  * Algorithm by Torben Mogensen, implementation by N. Devillard.
126  * This code in public domain.
127  * Adapted to SHOGUN by Heiko Strathmann
128  */
129  int32_t i;
130  int32_t less;
131  int32_t greater;
132  int32_t equal;
133  float64_t min;
134  float64_t max;
135  float64_t guess;
136  float64_t maxltguess;
137  float64_t mingtguess;
138  min=max=values[0];
139  for (i=1; i<values.vlen; i++)
140  {
141  if (values[i]<min)
142  min=values[i];
143  if (values[i]>max)
144  max=values[i];
145  }
146  while (1)
147  {
148  guess=(min+max)/2;
149  less=0;
150  greater=0;
151  equal=0;
152  maxltguess=min;
153  mingtguess=max;
154  for (i=0; i<values.vlen; i++)
155  {
156  if (values[i]<guess)
157  {
158  less++;
159  if (values[i]>maxltguess)
160  maxltguess=values[i];
161  }
162  else if (values[i]>guess)
163  {
164  greater++;
165  if (values[i]<mingtguess)
166  mingtguess=values[i];
167  }
168  else
169  equal++;
170  }
171  if (less<=(values.vlen+1)/2&&greater<=(values.vlen+1)/2)
172  break;
173  else if (less>greater)
174  max=maxltguess;
175  else
176  min=mingtguess;
177  }
178 
179  if (less>=(values.vlen+1)/2)
180  result=maxltguess;
181  else if (less+equal>=(values.vlen+1)/2)
182  result=guess;
183  else
184  result=mingtguess;
185  }
186  else
187  {
188  /* copy vector and do recursive call which modifies copy */
189  SGVector<float64_t> copy(values.vlen);
190  memcpy(copy.vector, values.vector, sizeof(float64_t)*values.vlen);
191  result=median(copy, true);
192  }
193  }
194 
195  return result;
196 }
197 
198 float64_t CStatistics::matrix_median(SGMatrix<float64_t> values,
199  bool modify, bool in_place)
200 {
201  /* create a vector that uses the matrix data, dont do reference counting */
202  SGVector<float64_t> as_vector(values.matrix,
203  values.num_rows*values.num_cols, false);
204 
205  /* return vector median method */
206  return median(as_vector, modify, in_place);
207 }
208 
209 
210 float64_t CStatistics::variance(SGVector<float64_t> values)
211 {
212  ASSERT(values.vlen>1)
213  ASSERT(values.vector)
214 
215  float64_t mean=CStatistics::mean(values);
216 
217  float64_t sum_squared_diff=0;
218  for (index_t i=0; i<values.vlen; ++i)
219  sum_squared_diff+=CMath::pow(values.vector[i]-mean, 2);
220 
221  return sum_squared_diff/(values.vlen-1);
222 }
223 
224 SGVector<float64_t> CStatistics::matrix_mean(SGMatrix<float64_t> values,
225  bool col_wise)
226 {
227  ASSERT(values.num_rows>0)
228  ASSERT(values.num_cols>0)
229  ASSERT(values.matrix)
230 
231  SGVector<float64_t> result;
232 
233  if (col_wise)
234  {
235  result=SGVector<float64_t>(values.num_cols);
236  for (index_t j=0; j<values.num_cols; ++j)
237  {
238  result[j]=0;
239  for (index_t i=0; i<values.num_rows; ++i)
240  result[j]+=values(i,j);
241 
242  result[j]/=values.num_rows;
243  }
244  }
245  else
246  {
247  result=SGVector<float64_t>(values.num_rows);
248  for (index_t j=0; j<values.num_rows; ++j)
249  {
250  result[j]=0;
251  for (index_t i=0; i<values.num_cols; ++i)
252  result[j]+=values(j,i);
253 
254  result[j]/=values.num_cols;
255  }
256  }
257 
258  return result;
259 }
260 
261 SGVector<float64_t> CStatistics::matrix_variance(SGMatrix<float64_t> values,
262  bool col_wise)
263 {
264  ASSERT(values.num_rows>0)
265  ASSERT(values.num_cols>0)
266  ASSERT(values.matrix)
267 
268  /* first compute mean */
269  SGVector<float64_t> mean=CStatistics::matrix_mean(values, col_wise);
270 
271  SGVector<float64_t> result;
272 
273  if (col_wise)
274  {
275  result=SGVector<float64_t>(values.num_cols);
276  for (index_t j=0; j<values.num_cols; ++j)
277  {
278  result[j]=0;
279  for (index_t i=0; i<values.num_rows; ++i)
280  result[j]+=CMath::pow(values(i,j)-mean[j], 2);
281 
282  result[j]/=(values.num_rows-1);
283  }
284  }
285  else
286  {
287  result=SGVector<float64_t>(values.num_rows);
288  for (index_t j=0; j<values.num_rows; ++j)
289  {
290  result[j]=0;
291  for (index_t i=0; i<values.num_cols; ++i)
292  result[j]+=CMath::pow(values(j,i)-mean[j], 2);
293 
294  result[j]/=(values.num_cols-1);
295  }
296  }
297 
298  return result;
299 }
300 
301 float64_t CStatistics::std_deviation(SGVector<float64_t> values)
302 {
303  return CMath::sqrt(variance(values));
304 }
305 
306 SGVector<float64_t> CStatistics::matrix_std_deviation(
307  SGMatrix<float64_t> values, bool col_wise)
308 {
309  SGVector<float64_t> var=CStatistics::matrix_variance(values, col_wise);
310  for (index_t i=0; i<var.vlen; ++i)
311  var[i]=CMath::sqrt(var[i]);
312 
313  return var;
314 }
315 
316 #ifdef HAVE_LAPACK
317 SGMatrix<float64_t> CStatistics::covariance_matrix(
318  SGMatrix<float64_t> observations, bool in_place)
319 {
320  SGMatrix<float64_t> centered=
321  in_place ?
322  observations :
323  SGMatrix<float64_t>(observations.num_rows,
324  observations.num_cols);
325 
326  if (!in_place)
327  {
328  memcpy(centered.matrix, observations.matrix,
329  sizeof(float64_t)*observations.num_rows*observations.num_cols);
330  }
331  centered.remove_column_mean();
332 
333  /* compute 1/(m-1) * X' * X */
334  SGMatrix<float64_t> cov=SGMatrix<float64_t>::matrix_multiply(centered,
335  centered, true, false, 1.0/(observations.num_rows-1));
336 
337  return cov;
338 }
339 #endif //HAVE_LAPACK
340 
341 float64_t CStatistics::confidence_intervals_mean(SGVector<float64_t> values,
342  float64_t alpha, float64_t& conf_int_low, float64_t& conf_int_up)
343 {
344  ASSERT(values.vlen>1)
345  ASSERT(values.vector)
346 
347  /* using one sided student t distribution evaluation */
348  alpha=alpha/2;
349 
350  /* degrees of freedom */
351  int32_t deg=values.vlen-1;
352 
353  /* compute absolute value of t-value */
354  float64_t t=CMath::abs(inverse_student_t(deg, alpha));
355 
356  /* values for calculating confidence interval */
357  float64_t std_dev=CStatistics::std_deviation(values);
358  float64_t mean=CStatistics::mean(values);
359 
360  /* compute confidence interval */
361  float64_t interval=t*std_dev/CMath::sqrt((float64_t)values.vlen);
362  conf_int_low=mean-interval;
363  conf_int_up=mean+interval;
364 
365  return mean;
366 }
367 
368 float64_t CStatistics::inverse_student_t(int32_t k, float64_t p)
369 {
370  float64_t t;
371  float64_t rk;
372  float64_t z;
373  int32_t rflg;
374  float64_t result;
375 
376  if (!(k>0 && greater(p, 0)) && less(p, 1))
377  {
378  SG_SERROR("CStatistics::inverse_student_t_distribution(): "
379  "Domain error\n");
380  }
381  rk=k;
382  if (greater(p, 0.25) && less(p, 0.75))
383  {
384  if (equal(p, 0.5))
385  {
386  result=0;
387  return result;
388  }
389  z=1.0-2.0*p;
390  z=inverse_incomplete_beta(0.5, 0.5*rk, CMath::abs(z));
391  t=CMath::sqrt(rk*z/(1.0-z));
392  if (less(p, 0.5))
393  {
394  t=-t;
395  }
396  result=t;
397  return result;
398  }
399  rflg=-1;
400  if (greater_equal(p, 0.5))
401  {
402  p=1.0-p;
403  rflg=1;
404  }
405  z=inverse_incomplete_beta(0.5*rk, 0.5, 2.0*p);
406  if (less(CMath::MAX_REAL_NUMBER*z, rk))
407  {
408  result=rflg*CMath::MAX_REAL_NUMBER;
409  return result;
410  }
411  t=CMath::sqrt(rk/z-rk);
412  result=rflg*t;
413  return result;
414 }
415 
416 float64_t CStatistics::inverse_incomplete_beta(float64_t a, float64_t b,
417  float64_t y)
418 {
419  float64_t aaa;
420  float64_t bbb;
421  float64_t y0;
422  float64_t d;
423  float64_t yyy;
424  float64_t x;
425  float64_t x0;
426  float64_t x1;
427  float64_t lgm;
428  float64_t yp;
429  float64_t di;
430  float64_t dithresh;
431  float64_t yl;
432  float64_t yh;
433  float64_t xt;
434  int32_t i;
435  int32_t rflg;
436  int32_t dir;
437  int32_t nflg;
438  int32_t mainlooppos;
439  int32_t ihalve;
440  int32_t ihalvecycle;
441  int32_t newt;
442  int32_t newtcycle;
443  int32_t breaknewtcycle;
444  int32_t breakihalvecycle;
445  float64_t result;
446 
447  i=0;
448  if (!(greater_equal(y, 0) && less_equal(y, 1)))
449  {
450  SG_SERROR("CStatistics::inverse_incomplete_beta(): "
451  "Domain error\n");
452  }
453 
454  /*
455  * special cases
456  */
457  if (equal(y, 0))
458  {
459  result=0;
460  return result;
461  }
462  if (equal(y, 1.0))
463  {
464  result=1;
465  return result;
466  }
467 
468  /*
469  * these initializations are not really necessary,
470  * but without them compiler complains about 'possibly uninitialized variables'.
471  */
472  dithresh=0;
473  rflg=0;
474  aaa=0;
475  bbb=0;
476  y0=0;
477  x=0;
478  yyy=0;
479  lgm=0;
480  dir=0;
481  di=0;
482 
483  /*
484  * normal initializations
485  */
486  x0=0.0;
487  yl=0.0;
488  x1=1.0;
489  yh=1.0;
490  nflg=0;
491  mainlooppos=0;
492  ihalve=1;
493  ihalvecycle=2;
494  newt=3;
495  newtcycle=4;
496  breaknewtcycle=5;
497  breakihalvecycle=6;
498 
499  /*
500  * main loop
501  */
502  for (;;)
503  {
504 
505  /*
506  * start
507  */
508  if (mainlooppos==0)
509  {
510  if (less_equal(a, 1.0) || less_equal(b, 1.0))
511  {
512  dithresh=1.0e-6;
513  rflg=0;
514  aaa=a;
515  bbb=b;
516  y0=y;
517  x=aaa/(aaa+bbb);
518  yyy=incomplete_beta(aaa, bbb, x);
519  mainlooppos=ihalve;
520  continue;
521  }
522  else
523  {
524  dithresh=1.0e-4;
525  }
526  yp=-inverse_normal_cdf(y);
527  if (greater(y, 0.5))
528  {
529  rflg=1;
530  aaa=b;
531  bbb=a;
532  y0=1.0-y;
533  yp=-yp;
534  }
535  else
536  {
537  rflg=0;
538  aaa=a;
539  bbb=b;
540  y0=y;
541  }
542  lgm=(yp*yp-3.0)/6.0;
543  x=2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
544  d=yp*CMath::sqrt(x+lgm)/x
545  -(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))
546  *(lgm+5.0/6.0-2.0/(3.0*x));
547  d=2.0*d;
548  if (less(d, CMath::log(CMath::MIN_REAL_NUMBER)))
549  {
550  x=0;
551  break;
552  }
553  x=aaa/(aaa+bbb*CMath::exp(d));
554  yyy=incomplete_beta(aaa, bbb, x);
555  yp=(yyy-y0)/y0;
556  if (less(CMath::abs(yp), 0.2))
557  {
558  mainlooppos=newt;
559  continue;
560  }
561  mainlooppos=ihalve;
562  continue;
563  }
564 
565  /*
566  * ihalve
567  */
568  if (mainlooppos==ihalve)
569  {
570  dir=0;
571  di=0.5;
572  i=0;
573  mainlooppos=ihalvecycle;
574  continue;
575  }
576 
577  /*
578  * ihalvecycle
579  */
580  if (mainlooppos==ihalvecycle)
581  {
582  if (i<=99)
583  {
584  if (i!=0)
585  {
586  x=x0+di*(x1-x0);
587  if (equal(x, 1.0))
588  {
589  x=1.0-CMath::MACHINE_EPSILON;
590  }
591  if (equal(x, 0.0))
592  {
593  di=0.5;
594  x=x0+di*(x1-x0);
595  if (equal(x, 0.0))
596  {
597  break;
598  }
599  }
600  yyy=incomplete_beta(aaa, bbb, x);
601  yp=(x1-x0)/(x1+x0);
602  if (less(CMath::abs(yp), dithresh))
603  {
604  mainlooppos=newt;
605  continue;
606  }
607  yp=(yyy-y0)/y0;
608  if (less(CMath::abs(yp), dithresh))
609  {
610  mainlooppos=newt;
611  continue;
612  }
613  }
614  if (less(yyy, y0))
615  {
616  x0=x;
617  yl=yyy;
618  if (dir<0)
619  {
620  dir=0;
621  di=0.5;
622  }
623  else
624  {
625  if (dir>3)
626  {
627  di=1.0-(1.0-di)*(1.0-di);
628  }
629  else
630  {
631  if (dir>1)
632  {
633  di=0.5*di+0.5;
634  }
635  else
636  {
637  di=(y0-yyy)/(yh-yl);
638  }
639  }
640  }
641  dir=dir+1;
642  if (greater(x0, 0.75))
643  {
644  if (rflg==1)
645  {
646  rflg=0;
647  aaa=a;
648  bbb=b;
649  y0=y;
650  }
651  else
652  {
653  rflg=1;
654  aaa=b;
655  bbb=a;
656  y0=1.0-y;
657  }
658  x=1.0-x;
659  yyy=incomplete_beta(aaa, bbb, x);
660  x0=0.0;
661  yl=0.0;
662  x1=1.0;
663  yh=1.0;
664  mainlooppos=ihalve;
665  continue;
666  }
667  }
668  else
669  {
670  x1=x;
671  if (rflg==1 && less(x1, CMath::MACHINE_EPSILON))
672  {
673  x=0.0;
674  break;
675  }
676  yh=yyy;
677  if (dir>0)
678  {
679  dir=0;
680  di=0.5;
681  }
682  else
683  {
684  if (dir<-3)
685  {
686  di=di*di;
687  }
688  else
689  {
690  if (dir<-1)
691  {
692  di=0.5*di;
693  }
694  else
695  {
696  di=(yyy-y0)/(yh-yl);
697  }
698  }
699  }
700  dir=dir-1;
701  }
702  i=i+1;
703  mainlooppos=ihalvecycle;
704  continue;
705  }
706  else
707  {
708  mainlooppos=breakihalvecycle;
709  continue;
710  }
711  }
712 
713  /*
714  * breakihalvecycle
715  */
716  if (mainlooppos==breakihalvecycle)
717  {
718  if (greater_equal(x0, 1.0))
719  {
720  x=1.0-CMath::MACHINE_EPSILON;
721  break;
722  }
723  if (less_equal(x, 0.0))
724  {
725  x=0.0;
726  break;
727  }
728  mainlooppos=newt;
729  continue;
730  }
731 
732  /*
733  * newt
734  */
735  if (mainlooppos==newt)
736  {
737  if (nflg!=0)
738  {
739  break;
740  }
741  nflg=1;
742  lgm=lgamma(aaa+bbb)-lgamma(aaa)-lgamma(bbb);
743  i=0;
744  mainlooppos=newtcycle;
745  continue;
746  }
747 
748  /*
749  * newtcycle
750  */
751  if (mainlooppos==newtcycle)
752  {
753  if (i<=7)
754  {
755  if (i!=0)
756  {
757  yyy=incomplete_beta(aaa, bbb, x);
758  }
759  if (less(yyy, yl))
760  {
761  x=x0;
762  yyy=yl;
763  }
764  else
765  {
766  if (greater(yyy, yh))
767  {
768  x=x1;
769  yyy=yh;
770  }
771  else
772  {
773  if (less(yyy, y0))
774  {
775  x0=x;
776  yl=yyy;
777  }
778  else
779  {
780  x1=x;
781  yh=yyy;
782  }
783  }
784  }
785  if (equal(x, 1.0) || equal(x, 0.0))
786  {
787  mainlooppos=breaknewtcycle;
788  continue;
789  }
790  d=(aaa-1.0)*CMath::log(x)+(bbb-1.0)*CMath::log(1.0-x)+lgm;
791  if (less(d, CMath::log(CMath::MIN_REAL_NUMBER)))
792  {
793  break;
794  }
795  if (greater(d, CMath::log(CMath::MAX_REAL_NUMBER)))
796  {
797  mainlooppos=breaknewtcycle;
798  continue;
799  }
800  d=CMath::exp(d);
801  d=(yyy-y0)/d;
802  xt=x-d;
803  if (less_equal(xt, x0))
804  {
805  yyy=(x-x0)/(x1-x0);
806  xt=x0+0.5*yyy*(x-x0);
807  if (less_equal(xt, 0.0))
808  {
809  mainlooppos=breaknewtcycle;
810  continue;
811  }
812  }
813  if (greater_equal(xt, x1))
814  {
815  yyy=(x1-x)/(x1-x0);
816  xt=x1-0.5*yyy*(x1-x);
817  if (greater_equal(xt, 1.0))
818  {
819  mainlooppos=breaknewtcycle;
820  continue;
821  }
822  }
823  x=xt;
824  if (less(CMath::abs(d/x), 128.0*CMath::MACHINE_EPSILON))
825  {
826  break;
827  }
828  i=i+1;
829  mainlooppos=newtcycle;
830  continue;
831  }
832  else
833  {
834  mainlooppos=breaknewtcycle;
835  continue;
836  }
837  }
838 
839  /*
840  * breaknewtcycle
841  */
842  if (mainlooppos==breaknewtcycle)
843  {
844  dithresh=256.0*CMath::MACHINE_EPSILON;
845  mainlooppos=ihalve;
846  continue;
847  }
848  }
849 
850  /*
851  * done
852  */
853  if (rflg!=0)
854  {
855  if (less_equal(x, CMath::MACHINE_EPSILON))
856  {
857  x=1.0-CMath::MACHINE_EPSILON;
858  }
859  else
860  {
861  x=1.0-x;
862  }
863  }
864  result=x;
865  return result;
866 }
867 
868 float64_t CStatistics::incomplete_beta(float64_t a, float64_t b, float64_t x)
869 {
870  float64_t t;
871  float64_t xc;
872  float64_t w;
873  float64_t y;
874  int32_t flag;
875  float64_t big;
876  float64_t biginv;
877  float64_t maxgam;
878  float64_t minlog;
879  float64_t maxlog;
880  float64_t result;
881 
882  big=4.503599627370496e15;
883  biginv=2.22044604925031308085e-16;
884  maxgam=171.624376956302725;
885  minlog=CMath::log(CMath::MIN_REAL_NUMBER);
886  maxlog=CMath::log(CMath::MAX_REAL_NUMBER);
887 
888  if (!(greater(a, 0) && greater(b, 0)))
889  {
890  SG_SERROR("CStatistics::incomplete_beta(): "
891  "Domain error\n");
892  }
893 
894  if (!(greater_equal(x, 0) && less_equal(x, 1)))
895  {
896  SG_SERROR("CStatistics::incomplete_beta(): "
897  "Domain error\n");
898  }
899 
900  if (equal(x, 0))
901  {
902  result=0;
903  return result;
904  }
905  if (equal(x, 1))
906  {
907  result=1;
908  return result;
909  }
910  flag=0;
911  if (less_equal(b*x, 1.0) && less_equal(x, 0.95))
912  {
913  result=ibetaf_incompletebetaps(a, b, x, maxgam);
914  return result;
915  }
916  w=1.0-x;
917  if (greater(x, a/(a+b)))
918  {
919  flag=1;
920  t=a;
921  a=b;
922  b=t;
923  xc=x;
924  x=w;
925  }
926  else
927  {
928  xc=w;
929  }
930  if ((flag==1 && less_equal(b*x, 1.0)) && less_equal(x, 0.95))
931  {
932  t=ibetaf_incompletebetaps(a, b, x, maxgam);
933  if (less_equal(t, CMath::MACHINE_EPSILON))
934  {
935  result=1.0-CMath::MACHINE_EPSILON;
936  }
937  else
938  {
939  result=1.0-t;
940  }
941  return result;
942  }
943  y=x*(a+b-2.0)-(a-1.0);
944  if (less(y, 0.0))
945  {
946  w=ibetaf_incompletebetafe(a, b, x, big, biginv);
947  }
948  else
949  {
950  w=ibetaf_incompletebetafe2(a, b, x, big, biginv)/xc;
951  }
952  y=a*CMath::log(x);
953  t=b*CMath::log(xc);
954  if ((less(a+b, maxgam) && less(CMath::abs(y), maxlog))
955  && less(CMath::abs(t), maxlog))
956  {
957  t=CMath::pow(xc, b);
958  t=t*CMath::pow(x, a);
959  t=t/a;
960  t=t*w;
961  t=t*(tgamma(a+b)/(tgamma(a)*tgamma(b)));
962  if (flag==1)
963  {
964  if (less_equal(t, CMath::MACHINE_EPSILON))
965  {
966  result=1.0-CMath::MACHINE_EPSILON;
967  }
968  else
969  {
970  result=1.0-t;
971  }
972  }
973  else
974  {
975  result=t;
976  }
977  return result;
978  }
979  y=y+t+lgamma(a+b)-lgamma(a)-lgamma(b);
980  y=y+CMath::log(w/a);
981  if (less(y, minlog))
982  {
983  t=0.0;
984  }
985  else
986  {
987  t=CMath::exp(y);
988  }
989  if (flag==1)
990  {
991  if (less_equal(t, CMath::MACHINE_EPSILON))
992  {
993  t=1.0-CMath::MACHINE_EPSILON;
994  }
995  else
996  {
997  t=1.0-t;
998  }
999  }
1000  result=t;
1001  return result;
1002 }
1003 
1004 float64_t CStatistics::inverse_normal_cdf(float64_t y0, float64_t mean,
1005  float64_t std_dev)
1006 {
1007  return inverse_normal_cdf(y0)*std_dev+mean;
1008 }
1009 
1010 float64_t CStatistics::inverse_normal_cdf(float64_t y0)
1011 {
1012  float64_t expm2;
1013  float64_t s2pi;
1014  float64_t x;
1015  float64_t y;
1016  float64_t z;
1017  float64_t y2;
1018  float64_t x0;
1019  float64_t x1;
1020  int32_t code;
1021  float64_t p0;
1022  float64_t q0;
1023  float64_t p1;
1024  float64_t q1;
1025  float64_t p2;
1026  float64_t q2;
1027  float64_t result;
1028 
1029  expm2=0.13533528323661269189;
1030  s2pi=2.50662827463100050242;
1031  if (less_equal(y0, 0))
1032  {
1033  result=-CMath::MAX_REAL_NUMBER;
1034  return result;
1035  }
1036  if (greater_equal(y0, 1))
1037  {
1038  result=CMath::MAX_REAL_NUMBER;
1039  return result;
1040  }
1041  code=1;
1042  y=y0;
1043  if (greater(y, 1.0-expm2))
1044  {
1045  y=1.0-y;
1046  code=0;
1047  }
1048  if (greater(y, expm2))
1049  {
1050  y=y-0.5;
1051  y2=y*y;
1052  p0=-59.9633501014107895267;
1053  p0=98.0010754185999661536+y2*p0;
1054  p0=-56.6762857469070293439+y2*p0;
1055  p0=13.9312609387279679503+y2*p0;
1056  p0=-1.23916583867381258016+y2*p0;
1057  q0=1;
1058  q0=1.95448858338141759834+y2*q0;
1059  q0=4.67627912898881538453+y2*q0;
1060  q0=86.3602421390890590575+y2*q0;
1061  q0=-225.462687854119370527+y2*q0;
1062  q0=200.260212380060660359+y2*q0;
1063  q0=-82.0372256168333339912+y2*q0;
1064  q0=15.9056225126211695515+y2*q0;
1065  q0=-1.18331621121330003142+y2*q0;
1066  x=y+y*y2*p0/q0;
1067  x=x*s2pi;
1068  result=x;
1069  return result;
1070  }
1071  x=CMath::sqrt(-2.0*CMath::log(y));
1072  x0=x-CMath::log(x)/x;
1073  z=1.0/x;
1074  if (less(x, 8.0))
1075  {
1076  p1=4.05544892305962419923;
1077  p1=31.5251094599893866154+z*p1;
1078  p1=57.1628192246421288162+z*p1;
1079  p1=44.0805073893200834700+z*p1;
1080  p1=14.6849561928858024014+z*p1;
1081  p1=2.18663306850790267539+z*p1;
1082  p1=-1.40256079171354495875*0.1+z*p1;
1083  p1=-3.50424626827848203418*0.01+z*p1;
1084  p1=-8.57456785154685413611*0.0001+z*p1;
1085  q1=1;
1086  q1=15.7799883256466749731+z*q1;
1087  q1=45.3907635128879210584+z*q1;
1088  q1=41.3172038254672030440+z*q1;
1089  q1=15.0425385692907503408+z*q1;
1090  q1=2.50464946208309415979+z*q1;
1091  q1=-1.42182922854787788574*0.1+z*q1;
1092  q1=-3.80806407691578277194*0.01+z*q1;
1093  q1=-9.33259480895457427372*0.0001+z*q1;
1094  x1=z*p1/q1;
1095  }
1096  else
1097  {
1098  p2=3.23774891776946035970;
1099  p2=6.91522889068984211695+z*p2;
1100  p2=3.93881025292474443415+z*p2;
1101  p2=1.33303460815807542389+z*p2;
1102  p2=2.01485389549179081538*0.1+z*p2;
1103  p2=1.23716634817820021358*0.01+z*p2;
1104  p2=3.01581553508235416007*0.0001+z*p2;
1105  p2=2.65806974686737550832*0.000001+z*p2;
1106  p2=6.23974539184983293730*0.000000001+z*p2;
1107  q2=1;
1108  q2=6.02427039364742014255+z*q2;
1109  q2=3.67983563856160859403+z*q2;
1110  q2=1.37702099489081330271+z*q2;
1111  q2=2.16236993594496635890*0.1+z*q2;
1112  q2=1.34204006088543189037*0.01+z*q2;
1113  q2=3.28014464682127739104*0.0001+z*q2;
1114  q2=2.89247864745380683936*0.000001+z*q2;
1115  q2=6.79019408009981274425*0.000000001+z*q2;
1116  x1=z*p2/q2;
1117  }
1118  x=x0-x1;
1119  if (code!=0)
1120  {
1121  x=-x;
1122  }
1123  result=x;
1124  return result;
1125 }
1126 
1127 float64_t CStatistics::ibetaf_incompletebetaps(float64_t a, float64_t b,
1128  float64_t x, float64_t maxgam)
1129 {
1130  float64_t s;
1131  float64_t t;
1132  float64_t u;
1133  float64_t v;
1134  float64_t n;
1135  float64_t t1;
1136  float64_t z;
1137  float64_t ai;
1138  float64_t result;
1139 
1140  ai=1.0/a;
1141  u=(1.0-b)*x;
1142  v=u/(a+1.0);
1143  t1=v;
1144  t=u;
1145  n=2.0;
1146  s=0.0;
1147  z=CMath::MACHINE_EPSILON*ai;
1148  while (greater(CMath::abs(v), z))
1149  {
1150  u=(n-b)*x/n;
1151  t=t*u;
1152  v=t/(a+n);
1153  s=s+v;
1154  n=n+1.0;
1155  }
1156  s=s+t1;
1157  s=s+ai;
1158  u=a*CMath::log(x);
1159  if (less(a+b, maxgam)
1160  && less(CMath::abs(u), CMath::log(CMath::MAX_REAL_NUMBER)))
1161  {
1162  t=tgamma(a+b)/(tgamma(a)*tgamma(b));
1163  s=s*t*CMath::pow(x, a);
1164  }
1165  else
1166  {
1167  t=lgamma(a+b)-lgamma(a)-lgamma(b)+u+CMath::log(s);
1168  if (less(t, CMath::log(CMath::MIN_REAL_NUMBER)))
1169  {
1170  s=0.0;
1171  }
1172  else
1173  {
1174  s=CMath::exp(t);
1175  }
1176  }
1177  result=s;
1178  return result;
1179 }
1180 
1181 float64_t CStatistics::ibetaf_incompletebetafe(float64_t a, float64_t b,
1182  float64_t x, float64_t big, float64_t biginv)
1183 {
1184  float64_t xk;
1185  float64_t pk;
1186  float64_t pkm1;
1187  float64_t pkm2;
1188  float64_t qk;
1189  float64_t qkm1;
1190  float64_t qkm2;
1191  float64_t k1;
1192  float64_t k2;
1193  float64_t k3;
1194  float64_t k4;
1195  float64_t k5;
1196  float64_t k6;
1197  float64_t k7;
1198  float64_t k8;
1199  float64_t r;
1200  float64_t t;
1201  float64_t ans;
1202  float64_t thresh;
1203  int32_t n;
1204  float64_t result;
1205 
1206  k1=a;
1207  k2=a+b;
1208  k3=a;
1209  k4=a+1.0;
1210  k5=1.0;
1211  k6=b-1.0;
1212  k7=k4;
1213  k8=a+2.0;
1214  pkm2=0.0;
1215  qkm2=1.0;
1216  pkm1=1.0;
1217  qkm1=1.0;
1218  ans=1.0;
1219  r=1.0;
1220  n=0;
1221  thresh=3.0*CMath::MACHINE_EPSILON;
1222  do
1223  {
1224  xk=-x*k1*k2/(k3*k4);
1225  pk=pkm1+pkm2*xk;
1226  qk=qkm1+qkm2*xk;
1227  pkm2=pkm1;
1228  pkm1=pk;
1229  qkm2=qkm1;
1230  qkm1=qk;
1231  xk=x*k5*k6/(k7*k8);
1232  pk=pkm1+pkm2*xk;
1233  qk=qkm1+qkm2*xk;
1234  pkm2=pkm1;
1235  pkm1=pk;
1236  qkm2=qkm1;
1237  qkm1=qk;
1238  if (not_equal(qk, 0))
1239  {
1240  r=pk/qk;
1241  }
1242  if (not_equal(r, 0))
1243  {
1244  t=CMath::abs((ans-r)/r);
1245  ans=r;
1246  }
1247  else
1248  {
1249  t=1.0;
1250  }
1251  if (less(t, thresh))
1252  {
1253  break;
1254  }
1255  k1=k1+1.0;
1256  k2=k2+1.0;
1257  k3=k3+2.0;
1258  k4=k4+2.0;
1259  k5=k5+1.0;
1260  k6=k6-1.0;
1261  k7=k7+2.0;
1262  k8=k8+2.0;
1263  if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1264  {
1265  pkm2=pkm2*biginv;
1266  pkm1=pkm1*biginv;
1267  qkm2=qkm2*biginv;
1268  qkm1=qkm1*biginv;
1269  }
1270  if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1271  {
1272  pkm2=pkm2*big;
1273  pkm1=pkm1*big;
1274  qkm2=qkm2*big;
1275  qkm1=qkm1*big;
1276  }
1277  n=n+1;
1278  }
1279  while (n!=300);
1280  result=ans;
1281  return result;
1282 }
1283 
1284 float64_t CStatistics::ibetaf_incompletebetafe2(float64_t a, float64_t b,
1285  float64_t x, float64_t big, float64_t biginv)
1286 {
1287  float64_t xk;
1288  float64_t pk;
1289  float64_t pkm1;
1290  float64_t pkm2;
1291  float64_t qk;
1292  float64_t qkm1;
1293  float64_t qkm2;
1294  float64_t k1;
1295  float64_t k2;
1296  float64_t k3;
1297  float64_t k4;
1298  float64_t k5;
1299  float64_t k6;
1300  float64_t k7;
1301  float64_t k8;
1302  float64_t r;
1303  float64_t t;
1304  float64_t ans;
1305  float64_t z;
1306  float64_t thresh;
1307  int32_t n;
1308  float64_t result;
1309 
1310  k1=a;
1311  k2=b-1.0;
1312  k3=a;
1313  k4=a+1.0;
1314  k5=1.0;
1315  k6=a+b;
1316  k7=a+1.0;
1317  k8=a+2.0;
1318  pkm2=0.0;
1319  qkm2=1.0;
1320  pkm1=1.0;
1321  qkm1=1.0;
1322  z=x/(1.0-x);
1323  ans=1.0;
1324  r=1.0;
1325  n=0;
1326  thresh=3.0*CMath::MACHINE_EPSILON;
1327  do
1328  {
1329  xk=-z*k1*k2/(k3*k4);
1330  pk=pkm1+pkm2*xk;
1331  qk=qkm1+qkm2*xk;
1332  pkm2=pkm1;
1333  pkm1=pk;
1334  qkm2=qkm1;
1335  qkm1=qk;
1336  xk=z*k5*k6/(k7*k8);
1337  pk=pkm1+pkm2*xk;
1338  qk=qkm1+qkm2*xk;
1339  pkm2=pkm1;
1340  pkm1=pk;
1341  qkm2=qkm1;
1342  qkm1=qk;
1343  if (not_equal(qk, 0))
1344  {
1345  r=pk/qk;
1346  }
1347  if (not_equal(r, 0))
1348  {
1349  t=CMath::abs((ans-r)/r);
1350  ans=r;
1351  }
1352  else
1353  {
1354  t=1.0;
1355  }
1356  if (less(t, thresh))
1357  {
1358  break;
1359  }
1360  k1=k1+1.0;
1361  k2=k2-1.0;
1362  k3=k3+2.0;
1363  k4=k4+2.0;
1364  k5=k5+1.0;
1365  k6=k6+1.0;
1366  k7=k7+2.0;
1367  k8=k8+2.0;
1368  if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1369  {
1370  pkm2=pkm2*biginv;
1371  pkm1=pkm1*biginv;
1372  qkm2=qkm2*biginv;
1373  qkm1=qkm1*biginv;
1374  }
1375  if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1376  {
1377  pkm2=pkm2*big;
1378  pkm1=pkm1*big;
1379  qkm2=qkm2*big;
1380  qkm1=qkm1*big;
1381  }
1382  n=n+1;
1383  }
1384  while (n!=300);
1385  result=ans;
1386  return result;
1387 }
1388 
1389 float64_t CStatistics::incomplete_gamma(float64_t a, float64_t x)
1390 {
1391  float64_t igammaepsilon;
1392  float64_t ans;
1393  float64_t ax;
1394  float64_t c;
1395  float64_t r;
1396  float64_t result;
1397 
1398  igammaepsilon=0.000000000000001;
1399  if (less_equal(x, 0) || less_equal(a, 0))
1400  {
1401  result=0;
1402  return result;
1403  }
1404  if (greater(x, 1) && greater(x, a))
1405  {
1406  result=1-incomplete_gamma_completed(a, x);
1407  return result;
1408  }
1409  ax=a*CMath::log(x)-x-lgamma(a);
1410  if (less(ax, -709.78271289338399))
1411  {
1412  result=0;
1413  return result;
1414  }
1415  ax=CMath::exp(ax);
1416  r=a;
1417  c=1;
1418  ans=1;
1419  do
1420  {
1421  r=r+1;
1422  c=c*x/r;
1423  ans=ans+c;
1424  }
1425  while (greater(c/ans, igammaepsilon));
1426  result=ans*ax/a;
1427  return result;
1428 }
1429 
1430 float64_t CStatistics::incomplete_gamma_completed(float64_t a, float64_t x)
1431 {
1432  float64_t igammaepsilon;
1433  float64_t igammabignumber;
1434  float64_t igammabignumberinv;
1435  float64_t ans;
1436  float64_t ax;
1437  float64_t c;
1438  float64_t yc;
1439  float64_t r;
1440  float64_t t;
1441  float64_t y;
1442  float64_t z;
1443  float64_t pk;
1444  float64_t pkm1;
1445  float64_t pkm2;
1446  float64_t qk;
1447  float64_t qkm1;
1448  float64_t qkm2;
1449  float64_t result;
1450 
1451  igammaepsilon=0.000000000000001;
1452  igammabignumber=4503599627370496.0;
1453  igammabignumberinv=2.22044604925031308085*0.0000000000000001;
1454  if (less_equal(x, 0) || less_equal(a, 0))
1455  {
1456  result=1;
1457  return result;
1458  }
1459  if (less(x, 1) || less(x, a))
1460  {
1461  result=1-incomplete_gamma(a, x);
1462  return result;
1463  }
1464  ax=a*CMath::log(x)-x-lgamma(a);
1465  if (less(ax, -709.78271289338399))
1466  {
1467  result=0;
1468  return result;
1469  }
1470  ax=CMath::exp(ax);
1471  y=1-a;
1472  z=x+y+1;
1473  c=0;
1474  pkm2=1;
1475  qkm2=x;
1476  pkm1=x+1;
1477  qkm1=z*x;
1478  ans=pkm1/qkm1;
1479  do
1480  {
1481  c=c+1;
1482  y=y+1;
1483  z=z+2;
1484  yc=y*c;
1485  pk=pkm1*z-pkm2*yc;
1486  qk=qkm1*z-qkm2*yc;
1487  if (not_equal(qk, 0))
1488  {
1489  r=pk/qk;
1490  t=CMath::abs((ans-r)/r);
1491  ans=r;
1492  }
1493  else
1494  {
1495  t=1;
1496  }
1497  pkm2=pkm1;
1498  pkm1=pk;
1499  qkm2=qkm1;
1500  qkm1=qk;
1501  if (greater(CMath::abs(pk), igammabignumber))
1502  {
1503  pkm2=pkm2*igammabignumberinv;
1504  pkm1=pkm1*igammabignumberinv;
1505  qkm2=qkm2*igammabignumberinv;
1506  qkm1=qkm1*igammabignumberinv;
1507  }
1508  }
1509  while (greater(t, igammaepsilon));
1510  result=ans*ax;
1511  return result;
1512 }
1513 
1514 float64_t CStatistics::gamma_cdf(float64_t x, float64_t a, float64_t b)
1515 {
1516  /* definition of wikipedia: incomplete gamma devised by true gamma */
1517  return incomplete_gamma(a, x/b);
1518 }
1519 
1520 float64_t CStatistics::inverse_gamma_cdf(float64_t p, float64_t a,
1521  float64_t b)
1522 {
1523  /* inverse of gamma(a,b) CDF is
1524  * inverse_incomplete_gamma_completed(a, 1. - p) * b */
1525  return inverse_incomplete_gamma_completed(a, 1-p)*b;
1526 }
1527 
1528 float64_t CStatistics::inverse_incomplete_gamma_completed(float64_t a,
1529  float64_t y0)
1530 {
1531  float64_t igammaepsilon;
1532  float64_t iinvgammabignumber;
1533  float64_t x0;
1534  float64_t x1;
1535  float64_t x;
1536  float64_t yl;
1537  float64_t yh;
1538  float64_t y;
1539  float64_t d;
1540  float64_t lgm;
1541  float64_t dithresh;
1542  int32_t i;
1543  int32_t dir;
1544  float64_t result;
1545 
1546  igammaepsilon=0.000000000000001;
1547  iinvgammabignumber=4503599627370496.0;
1548  x0=iinvgammabignumber;
1549  yl=0;
1550  x1=0;
1551  yh=1;
1552  dithresh=5*igammaepsilon;
1553  d=1/(9*a);
1554  y=1-d-inverse_normal_cdf(y0)*CMath::sqrt(d);
1555  x=a*y*y*y;
1556  lgm=lgamma(a);
1557  i=0;
1558  while (i<10)
1559  {
1560  if (greater(x, x0) || less(x, x1))
1561  {
1562  d=0.0625;
1563  break;
1564  }
1565  y=incomplete_gamma_completed(a, x);
1566  if (less(y, yl) || greater(y, yh))
1567  {
1568  d=0.0625;
1569  break;
1570  }
1571  if (less(y, y0))
1572  {
1573  x0=x;
1574  yl=y;
1575  }
1576  else
1577  {
1578  x1=x;
1579  yh=y;
1580  }
1581  d=(a-1)*CMath::log(x)-x-lgm;
1582  if (less(d, -709.78271289338399))
1583  {
1584  d=0.0625;
1585  break;
1586  }
1587  d=-CMath::exp(d);
1588  d=(y-y0)/d;
1589  if (less(CMath::abs(d/x), igammaepsilon))
1590  {
1591  result=x;
1592  return result;
1593  }
1594  x=x-d;
1595  i=i+1;
1596  }
1597  if (equal(x0, iinvgammabignumber))
1598  {
1599  if (less_equal(x, 0))
1600  {
1601  x=1;
1602  }
1603  while (equal(x0, iinvgammabignumber))
1604  {
1605  x=(1+d)*x;
1606  y=incomplete_gamma_completed(a, x);
1607  if (less(y, y0))
1608  {
1609  x0=x;
1610  yl=y;
1611  break;
1612  }
1613  d=d+d;
1614  }
1615  }
1616  d=0.5;
1617  dir=0;
1618  i=0;
1619  while (i<400)
1620  {
1621  x=x1+d*(x0-x1);
1622  y=incomplete_gamma_completed(a, x);
1623  lgm=(x0-x1)/(x1+x0);
1624  if (less(CMath::abs(lgm), dithresh))
1625  {
1626  break;
1627  }
1628  lgm=(y-y0)/y0;
1629  if (less(CMath::abs(lgm), dithresh))
1630  {
1631  break;
1632  }
1633  if (less_equal(x, 0.0))
1634  {
1635  break;
1636  }
1637  if (greater_equal(y, y0))
1638  {
1639  x1=x;
1640  yh=y;
1641  if (dir<0)
1642  {
1643  dir=0;
1644  d=0.5;
1645  }
1646  else
1647  {
1648  if (dir>1)
1649  {
1650  d=0.5*d+0.5;
1651  }
1652  else
1653  {
1654  d=(y0-yl)/(yh-yl);
1655  }
1656  }
1657  dir=dir+1;
1658  }
1659  else
1660  {
1661  x0=x;
1662  yl=y;
1663  if (dir>0)
1664  {
1665  dir=0;
1666  d=0.5;
1667  }
1668  else
1669  {
1670  if (dir<-1)
1671  {
1672  d=0.5*d;
1673  }
1674  else
1675  {
1676  d=(y0-yl)/(yh-yl);
1677  }
1678  }
1679  dir=dir-1;
1680  }
1681  i=i+1;
1682  }
1683  result=x;
1684  return result;
1685 }
1686 
1687 float64_t CStatistics::normal_cdf(float64_t x, float64_t std_dev)
1688 {
1689  return 0.5*(error_function(x/std_dev/1.41421356237309504880)+1);
1690 }
1691 
1692 float64_t CStatistics::lnormal_cdf(float64_t x)
1693 {
1694  float64_t result;
1695 
1696  if (x<-10.0)
1697  {
1698  float64_t x2=x*x;
1699  float64_t s=1.0-1.0/x2*(1.0-3.0/x2*(1.0-5.0/x2*(1.0-7.0/x2)));
1700  result=-0.5*CMath::log(2*CMath::PI)-x2*0.5-CMath::log(-x)+CMath::log(s);
1701  }
1702  else
1703  result=CMath::log(normal_cdf(x));
1704 
1705  return result;
1706 }
1707 
1708 float64_t CStatistics::error_function(float64_t x)
1709 {
1710  float64_t xsq;
1711  float64_t s;
1712  float64_t p;
1713  float64_t q;
1714  float64_t result;
1715 
1716  s=CMath::sign(x);
1717  x=CMath::abs(x);
1718  if (less(x, 0.5))
1719  {
1720  xsq=x*x;
1721  p=0.007547728033418631287834;
1722  p=0.288805137207594084924010+xsq*p;
1723  p=14.3383842191748205576712+xsq*p;
1724  p=38.0140318123903008244444+xsq*p;
1725  p=3017.82788536507577809226+xsq*p;
1726  p=7404.07142710151470082064+xsq*p;
1727  p=80437.3630960840172832162+xsq*p;
1728  q=0.0;
1729  q=1.00000000000000000000000+xsq*q;
1730  q=38.0190713951939403753468+xsq*q;
1731  q=658.070155459240506326937+xsq*q;
1732  q=6379.60017324428279487120+xsq*q;
1733  q=34216.5257924628539769006+xsq*q;
1734  q=80437.3630960840172826266+xsq*q;
1735  result=s*1.1283791670955125738961589031*x*p/q;
1736  return result;
1737  }
1738  if (greater_equal(x, 10))
1739  {
1740  result=s;
1741  return result;
1742  }
1743  result=s*(1-error_function_complement(x));
1744  return result;
1745 }
1746 
1747 float64_t CStatistics::error_function_complement(float64_t x)
1748 {
1749  float64_t p;
1750  float64_t q;
1751  float64_t result;
1752 
1753  if (less(x, 0))
1754  {
1755  result=2-error_function_complement(-x);
1756  return result;
1757  }
1758  if (less(x, 0.5))
1759  {
1760  result=1.0-error_function(x);
1761  return result;
1762  }
1763  if (greater_equal(x, 10))
1764  {
1765  result=0;
1766  return result;
1767  }
1768  p=0.0;
1769  p=0.5641877825507397413087057563+x*p;
1770  p=9.675807882987265400604202961+x*p;
1771  p=77.08161730368428609781633646+x*p;
1772  p=368.5196154710010637133875746+x*p;
1773  p=1143.262070703886173606073338+x*p;
1774  p=2320.439590251635247384768711+x*p;
1775  p=2898.0293292167655611275846+x*p;
1776  p=1826.3348842295112592168999+x*p;
1777  q=1.0;
1778  q=17.14980943627607849376131193+x*q;
1779  q=137.1255960500622202878443578+x*q;
1780  q=661.7361207107653469211984771+x*q;
1781  q=2094.384367789539593790281779+x*q;
1782  q=4429.612803883682726711528526+x*q;
1783  q=6089.5424232724435504633068+x*q;
1784  q=4958.82756472114071495438422+x*q;
1785  q=1826.3348842295112595576438+x*q;
1786  result=CMath::exp(-x*x)*p/q;
1787  return result;
1788 }
1789 
1790 SGVector<float64_t> CStatistics::fishers_exact_test_for_multiple_2x3_tables(
1791  SGMatrix<float64_t> tables)
1792 {
1793  SGMatrix<float64_t> table(NULL, 2, 3, false);
1794  int32_t len=tables.num_cols/3;
1795 
1796  SGVector<float64_t> v(len);
1797  for (int32_t i=0; i<len; i++)
1798  {
1799  table.matrix=&tables.matrix[2*3*i];
1800  v.vector[i]=fishers_exact_test_for_2x3_table(table);
1801  }
1802  return v;
1803 }
1804 
1805 float64_t CStatistics::fishers_exact_test_for_2x3_table(
1806  SGMatrix<float64_t> table)
1807 {
1808  ASSERT(table.num_rows==2)
1809  ASSERT(table.num_cols==3)
1810 
1811  int32_t m_len=3+2;
1812  float64_t* m=SG_MALLOC(float64_t, 3+2);
1813  m[0]=table.matrix[0]+table.matrix[2]+table.matrix[4];
1814  m[1]=table.matrix[1]+table.matrix[3]+table.matrix[5];
1815  m[2]=table.matrix[0]+table.matrix[1];
1816  m[3]=table.matrix[2]+table.matrix[3];
1817  m[4]=table.matrix[4]+table.matrix[5];
1818 
1819  float64_t n=SGVector<float64_t>::sum(m, m_len)/2.0;
1820  int32_t x_len=2*3*CMath::sq(SGVector<float64_t>::max(m, m_len));
1821  float64_t* x=SG_MALLOC(float64_t, x_len);
1822  SGVector<float64_t>::fill_vector(x, x_len, 0.0);
1823 
1824  float64_t log_nom=0.0;
1825  for (int32_t i=0; i<3+2; i++)
1826  log_nom+=lgamma(m[i]+1);
1827  log_nom-=lgamma(n+1.0);
1828 
1829  float64_t log_denomf=0;
1830  floatmax_t log_denom=0;
1831 
1832  for (int32_t i=0; i<3*2; i++)
1833  {
1834  log_denom+=lgammal((floatmax_t)table.matrix[i]+1);
1835  log_denomf+=lgammal((floatmax_t)table.matrix[i]+1);
1836  }
1837 
1838  floatmax_t prob_table_log=log_nom-log_denom;
1839 
1840  int32_t dim1=CMath::min(m[0], m[2]);
1841 
1842  //traverse all possible tables with given m
1843  int32_t counter=0;
1844  for (int32_t k=0; k<=dim1; k++)
1845  {
1846  for (int32_t l=CMath::max(0.0, m[0]-m[4]-k);
1847  l<=CMath::min(m[0]-k, m[3]); l++)
1848  {
1849  x[0+0*2+counter*2*3]=k;
1850  x[0+1*2+counter*2*3]=l;
1851  x[0+2*2+counter*2*3]=m[0]-x[0+0*2+counter*2*3]-x[0+1*2+counter*2*3];
1852  x[1+0*2+counter*2*3]=m[2]-x[0+0*2+counter*2*3];
1853  x[1+1*2+counter*2*3]=m[3]-x[0+1*2+counter*2*3];
1854  x[1+2*2+counter*2*3]=m[4]-x[0+2*2+counter*2*3];
1855 
1856  counter++;
1857  }
1858  }
1859 
1860 //#define DEBUG_FISHER_TABLE
1861 #ifdef DEBUG_FISHER_TABLE
1862  SG_SPRINT("counter=%d\n", counter)
1863  SG_SPRINT("dim1=%d\n", dim1)
1864  SG_SPRINT("l=%g...%g\n", CMath::max(0.0,m[0]-m[4]-0), CMath::min(m[0]-0, m[3]))
1865  SG_SPRINT("n=%g\n", n)
1866  SG_SPRINT("prob_table_log=%.18Lg\n", prob_table_log)
1867  SG_SPRINT("log_denomf=%.18g\n", log_denomf)
1868  SG_SPRINT("log_denom=%.18Lg\n", log_denom)
1869  SG_SPRINT("log_nom=%.18g\n", log_nom)
1870  display_vector(m, m_len, "marginals");
1871  display_vector(x, 2*3*counter, "x");
1872 #endif // DEBUG_FISHER_TABLE
1873 
1874  floatmax_t* log_denom_vec=SG_MALLOC(floatmax_t, counter);
1875  SGVector<floatmax_t>::fill_vector(log_denom_vec, counter, (floatmax_t)0.0);
1876 
1877  for (int32_t k=0; k<counter; k++)
1878  {
1879  for (int32_t j=0; j<3; j++)
1880  {
1881  for (int32_t i=0; i<2; i++)
1882  log_denom_vec[k]+=lgammal(x[i+j*2+k*2*3]+1.0);
1883  }
1884  }
1885 
1886  for (int32_t i=0; i<counter; i++)
1887  log_denom_vec[i]=log_nom-log_denom_vec[i];
1888 
1889 #ifdef DEBUG_FISHER_TABLE
1890  display_vector(log_denom_vec, counter, "log_denom_vec");
1891 #endif // DEBUG_FISHER_TABLE
1892 
1893  float64_t nonrand_p=-CMath::INFTY;
1894  for (int32_t i=0; i<counter; i++)
1895  {
1896  if (log_denom_vec[i]<=prob_table_log)
1897  nonrand_p=CMath::logarithmic_sum(nonrand_p, log_denom_vec[i]);
1898  }
1899 
1900 #ifdef DEBUG_FISHER_TABLE
1901  SG_SPRINT("nonrand_p=%.18g\n", nonrand_p)
1902  SG_SPRINT("exp_nonrand_p=%.18g\n", CMath::exp(nonrand_p))
1903 #endif // DEBUG_FISHER_TABLE
1904  nonrand_p=CMath::exp(nonrand_p);
1905 
1906  SG_FREE(log_denom_vec);
1907  SG_FREE(x);
1908  SG_FREE(m);
1909 
1910  return nonrand_p;
1911 }
1912 
1913 float64_t CStatistics::mutual_info(float64_t* p1, float64_t* p2, int32_t len)
1914 {
1915  double e=0;
1916 
1917  for (int32_t i=0; i<len; i++)
1918  for (int32_t j=0; j<len; j++)
1919  e+=exp(p2[j*len+i])*(p2[j*len+i]-p1[i]-p1[j]);
1920 
1921  return (float64_t)e;
1922 }
1923 
1924 float64_t CStatistics::relative_entropy(float64_t* p, float64_t* q, int32_t len)
1925 {
1926  double e=0;
1927 
1928  for (int32_t i=0; i<len; i++)
1929  e+=exp(p[i])*(p[i]-q[i]);
1930 
1931  return (float64_t)e;
1932 }
1933 
1934 float64_t CStatistics::entropy(float64_t* p, int32_t len)
1935 {
1936  double e=0;
1937 
1938  for (int32_t i=0; i<len; i++)
1939  e-=exp(p[i])*p[i];
1940 
1941  return (float64_t)e;
1942 }
1943 
1944 SGVector<int32_t> CStatistics::sample_indices(int32_t sample_size, int32_t N)
1945 {
1946  REQUIRE(sample_size<N,
1947  "sample size should be less than number of indices\n");
1948  int32_t* idxs=SG_MALLOC(int32_t,N);
1949  int32_t i, rnd;
1950  int32_t* permuted_idxs=SG_MALLOC(int32_t,sample_size);
1951 
1952  // reservoir sampling
1953  for (i=0; i<N; i++)
1954  idxs[i]=i;
1955  for (i=0; i<sample_size; i++)
1956  permuted_idxs[i]=idxs[i];
1957  for (i=sample_size; i<N; i++)
1958  {
1959  rnd=CMath::random(1, i);
1960  if (rnd<sample_size)
1961  permuted_idxs[rnd]=idxs[i];
1962  }
1963  SG_FREE(idxs);
1964 
1965  SGVector<int32_t> result=SGVector<int32_t>(permuted_idxs, sample_size);
1966  result.qsort();
1967  return result;
1968 }
1969 
1970 float64_t CStatistics::dlgamma(float64_t x)
1971 {
1972  float64_t result=0.0;
1973 
1974  if (x<0.0)
1975  {
1976  // use reflection formula
1977  x=1.0-x;
1978  result=CMath::PI/CMath::tan(CMath::PI*x);
1979  }
1980 
1981  // make x>7 for approximation
1982  // (use reccurent formula: psi(x+1) = psi(x) + 1/x)
1983  while (x<=7.0)
1984  {
1985  result-=1.0/x;
1986  x++;
1987  }
1988 
1989  // perform approximation
1990  x-=0.5;
1991  result+=log(x);
1992 
1993  float64_t coeff[10]={
1994  0.04166666666666666667,
1995  -0.00729166666666666667,
1996  0.00384424603174603175,
1997  -0.00413411458333333333,
1998  0.00756096117424242424,
1999  -0.02108249687595390720,
2000  0.08332316080729166666,
2001  -0.44324627670587277880,
2002  3.05393103044765369366,
2003  -26.45616165999210241989};
2004 
2005  float64_t power=1.0;
2006  float64_t ix2=1.0/CMath::sq(x);
2007 
2008  // perform approximation
2009  for (index_t i=0; i<10; i++)
2010  {
2011  power*=ix2;
2012  result+=coeff[i]*power;
2013  }
2014 
2015  return result;
2016 }
2017 
2018 #ifdef HAVE_EIGEN3
2019 float64_t CStatistics::log_det(SGMatrix<float64_t> m)
2020 {
2021  /* map the matrix to eigen3 to perform cholesky */
2022  Map<MatrixXd> M(m.matrix, m.num_rows, m.num_cols);
2023 
2024  /* computing the cholesky decomposition */
2025  LLT<MatrixXd> llt;
2026  llt.compute(M);
2027 
2028  /* the lower triangular matrix */
2029  MatrixXd l = llt.matrixL();
2030 
2031  /* calculate the log-determinant */
2032  VectorXd diag = l.diagonal();
2033  float64_t retval = 0.0;
2034  for( int32_t i = 0; i < diag.rows(); ++i ) {
2035  retval += log(diag(i));
2036  }
2037  retval *= 2;
2038 
2039  return retval;
2040 }
2041 
2042 float64_t CStatistics::log_det(const SGSparseMatrix<float64_t> m)
2043 {
2044  typedef SparseMatrix<float64_t> MatrixType;
2045  const MatrixType &M=EigenSparseUtil<float64_t>::toEigenSparse(m);
2046 
2047  SimplicialLLT<MatrixType> llt;
2048 
2049  // factorize using cholesky with amd permutation
2050  llt.compute(M);
2051  MatrixType L=llt.matrixL();
2052 
2053  // calculate the log-determinant
2054  float64_t retval=0.0;
2055  for( index_t i=0; i<M.rows(); ++i )
2056  retval+=log(L.coeff(i,i));
2057  retval*=2;
2058 
2059  return retval;
2060 }
2061 
2062 SGMatrix<float64_t> CStatistics::sample_from_gaussian(SGVector<float64_t> mean,
2063  SGMatrix<float64_t> cov, int32_t N, bool precision_matrix)
2064 {
2065  REQUIRE(cov.num_rows>0,
2066  "CStatistics::sample_from_gaussian(): \
2067  Number of covariance rows must be positive!\n");
2068  REQUIRE(cov.num_cols>0,
2069  "CStatistics::sample_from_gaussian(): \
2070  Number of covariance cols must be positive!\n");
2071  REQUIRE(cov.matrix,
2072  "CStatistics::sample_from_gaussian(): \
2073  Covariance is not initialized!\n");
2074  REQUIRE(cov.num_rows==cov.num_cols,
2075  "CStatistics::sample_from_gaussian(): \
2076  Covariance should be square matrix!\n");
2077  REQUIRE(mean.vlen==cov.num_rows,
2078  "CStatistics::sample_from_gaussian(): \
2079  Mean and covariance dimension mismatch!\n");
2080 
2081  int32_t dim=mean.vlen;
2082  Map<VectorXd> mu(mean.vector, mean.vlen);
2083  Map<MatrixXd> c(cov.matrix, cov.num_rows, cov.num_cols);
2084 
2085  // generate samples, z, from N(0, I), DxN
2086  SGMatrix<float64_t> S(dim, N);
2087  for( int32_t j=0; j<N; ++j )
2088  for( int32_t i=0; i<dim; ++i )
2089  S(i,j)=CMath::randn_double();
2090 
2091  // the cholesky factorization c=L*U
2092  MatrixXd U=c.llt().matrixU();
2093 
2094  // generate samples, x, from N(mean, cov) or N(mean, cov^-1)
2095  // return samples of dimension NxD
2096  if( precision_matrix )
2097  {
2098  // here we have U*x=z, to solve this, we use cholesky again
2099  Map<MatrixXd> s(S.matrix, S.num_rows, S.num_cols);
2100  LDLT<MatrixXd> ldlt;
2101  ldlt.compute(U);
2102  s=ldlt.solve(s);
2103  }
2104 
2105  SGMatrix<float64_t>::transpose_matrix(S.matrix, S.num_rows, S.num_cols);
2106 
2107  if( !precision_matrix )
2108  {
2109  // here we need to find x=L*z, so, x'=z'*L' i.e. x'=z'*U
2110  Map<MatrixXd> s(S.matrix, S.num_rows, S.num_cols);
2111  s=s*U;
2112  }
2113 
2114  // add the mean
2115  Map<MatrixXd> x(S.matrix, S.num_rows, S.num_cols);
2116  for( int32_t i=0; i<N; ++i )
2117  x.row(i)+=mu;
2118 
2119  return S;
2120 }
2121 
2122 SGMatrix<float64_t> CStatistics::sample_from_gaussian(SGVector<float64_t> mean,
2123  SGSparseMatrix<float64_t> cov, int32_t N, bool precision_matrix)
2124 {
2125  REQUIRE(cov.num_vectors>0,
2126  "CStatistics::sample_from_gaussian(): \
2127  Number of covariance rows must be positive!\n");
2128  REQUIRE(cov.num_features>0,
2129  "CStatistics::sample_from_gaussian(): \
2130  Number of covariance cols must be positive!\n");
2131  REQUIRE(cov.sparse_matrix,
2132  "CStatistics::sample_from_gaussian(): \
2133  Covariance is not initialized!\n");
2134  REQUIRE(cov.num_vectors==cov.num_features,
2135  "CStatistics::sample_from_gaussian(): \
2136  Covariance should be square matrix!\n");
2137  REQUIRE(mean.vlen==cov.num_vectors,
2138  "CStatistics::sample_from_gaussian(): \
2139  Mean and covariance dimension mismatch!\n");
2140 
2141  int32_t dim=mean.vlen;
2142  Map<VectorXd> mu(mean.vector, mean.vlen);
2143 
2144  typedef SparseMatrix<float64_t> MatrixType;
2145  const MatrixType &c=EigenSparseUtil<float64_t>::toEigenSparse(cov);
2146 
2147  SimplicialLLT<MatrixType> llt;
2148 
2149  // generate samples, z, from N(0, I), DxN
2150  SGMatrix<float64_t> S(dim, N);
2151  for( int32_t j=0; j<N; ++j )
2152  for( int32_t i=0; i<dim; ++i )
2153  S(i,j)=CMath::randn_double();
2154 
2155  Map<MatrixXd> s(S.matrix, S.num_rows, S.num_cols);
2156 
2157  // the cholesky factorization P*c*P^-1 = LP*UP, with LP=P*L, UP=U*P^-1
2158  llt.compute(c);
2159  MatrixType LP=llt.matrixL();
2160  MatrixType UP=llt.matrixU();
2161 
2162  // generate samples, x, from N(mean, cov) or N(mean, cov^-1)
2163  // return samples of dimension NxD
2164  if( precision_matrix )
2165  {
2166  // here we have UP*xP=z, to solve this, we use cholesky again
2167  SimplicialLLT<MatrixType> lltUP;
2168  lltUP.compute(UP);
2169  s=lltUP.solve(s);
2170  }
2171  else
2172  {
2173  // here we need to find xP=LP*z
2174  s=LP*s;
2175  }
2176 
2177  // permute the samples back with x=P^-1*xP
2178  s=llt.permutationPinv()*s;
2179 
2180  SGMatrix<float64_t>::transpose_matrix(S.matrix, S.num_rows, S.num_cols);
2181  // add the mean
2182  Map<MatrixXd> x(S.matrix, S.num_rows, S.num_cols);
2183  for( int32_t i=0; i<N; ++i )
2184  x.row(i)+=mu;
2185 
2186  return S;
2187 }
2188 
2189 #endif //HAVE_EIGEN3
2190 
2191 CStatistics::SigmoidParamters CStatistics::fit_sigmoid(SGVector<float64_t> scores)
2192 {
2193  SG_SDEBUG("entering CStatistics::fit_sigmoid()\n")
2194 
2195  REQUIRE(scores.vector, "CStatistics::fit_sigmoid() requires "
2196  "scores vector!\n");
2197 
2198  /* count prior0 and prior1 if needed */
2199  int32_t prior0=0;
2200  int32_t prior1=0;
2201  SG_SDEBUG("counting number of positive and negative labels\n")
2202  {
2203  for (index_t i=0; i<scores.vlen; ++i)
2204  {
2205  if (scores[i]>0)
2206  prior1++;
2207  else
2208  prior0++;
2209  }
2210  }
2211  SG_SDEBUG("%d pos; %d neg\n", prior1, prior0)
2212 
2213  /* parameter setting */
2214  /* maximum number of iterations */
2215  index_t maxiter=100;
2216 
2217  /* minimum step taken in line search */
2218  float64_t minstep=1E-10;
2219 
2220  /* for numerically strict pd of hessian */
2221  float64_t sigma=1E-12;
2222  float64_t eps=1E-5;
2223 
2224  /* construct target support */
2225  float64_t hiTarget=(prior1+1.0)/(prior1+2.0);
2226  float64_t loTarget=1/(prior0+2.0);
2227  index_t length=prior1+prior0;
2228 
2229  SGVector<float64_t> t(length);
2230  for (index_t i=0; i<length; ++i)
2231  {
2232  if (scores[i]>0)
2233  t[i]=hiTarget;
2234  else
2235  t[i]=loTarget;
2236  }
2237 
2238  /* initial Point and Initial Fun Value */
2239  /* result parameters of sigmoid */
2240  float64_t a=0;
2241  float64_t b=CMath::log((prior0+1.0)/(prior1+1.0));
2242  float64_t fval=0.0;
2243 
2244  for (index_t i=0; i<length; ++i)
2245  {
2246  float64_t fApB=scores[i]*a+b;
2247  if (fApB>=0)
2248  fval+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
2249  else
2250  fval+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
2251  }
2252 
2253  index_t it;
2254  float64_t g1;
2255  float64_t g2;
2256  for (it=0; it<maxiter; ++it)
2257  {
2258  SG_SDEBUG("Iteration %d, a=%f, b=%f, fval=%f\n", it, a, b, fval)
2259 
2260  /* Update Gradient and Hessian (use H' = H + sigma I) */
2261  float64_t h11=sigma; //Numerically ensures strict PD
2262  float64_t h22=h11;
2263  float64_t h21=0;
2264  g1=0;
2265  g2=0;
2266 
2267  for (index_t i=0; i<length; ++i)
2268  {
2269  float64_t fApB=scores[i]*a+b;
2270  float64_t p;
2271  float64_t q;
2272  if (fApB>=0)
2273  {
2274  p=CMath::exp(-fApB)/(1.0+CMath::exp(-fApB));
2275  q=1.0/(1.0+CMath::exp(-fApB));
2276  }
2277  else
2278  {
2279  p=1.0/(1.0+CMath::exp(fApB));
2280  q=CMath::exp(fApB)/(1.0+CMath::exp(fApB));
2281  }
2282 
2283  float64_t d2=p*q;
2284  h11+=scores[i]*scores[i]*d2;
2285  h22+=d2;
2286  h21+=scores[i]*d2;
2287  float64_t d1=t[i]-p;
2288  g1+=scores[i]*d1;
2289  g2+=d1;
2290  }
2291 
2292  /* Stopping Criteria */
2293  if (CMath::abs(g1)<eps && CMath::abs(g2)<eps)
2294  break;
2295 
2296  /* Finding Newton direction: -inv(H') * g */
2297  float64_t det=h11*h22-h21*h21;
2298  float64_t dA=-(h22*g1-h21*g2)/det;
2299  float64_t dB=-(-h21*g1+h11*g2)/det;
2300  float64_t gd=g1*dA+g2*dB;
2301 
2302  /* Line Search */
2303  float64_t stepsize=1;
2304 
2305  while (stepsize>=minstep)
2306  {
2307  float64_t newA=a+stepsize*dA;
2308  float64_t newB=b+stepsize*dB;
2309 
2310  /* New function value */
2311  float64_t newf=0.0;
2312  for (index_t i=0; i<length; ++i)
2313  {
2314  float64_t fApB=scores[i]*newA+newB;
2315  if (fApB>=0)
2316  newf+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
2317  else
2318  newf+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
2319  }
2320 
2321  /* Check sufficient decrease */
2322  if (newf<fval+0.0001*stepsize*gd)
2323  {
2324  a=newA;
2325  b=newB;
2326  fval=newf;
2327  break;
2328  }
2329  else
2330  stepsize=stepsize/2.0;
2331  }
2332 
2333  if (stepsize<minstep)
2334  {
2335  SG_SWARNING("CStatistics::fit_sigmoid(): line search fails, A=%f, "
2336  "B=%f, g1=%f, g2=%f, dA=%f, dB=%f, gd=%f\n",
2337  a, b, g1, g2, dA, dB, gd);
2338  }
2339  }
2340 
2341  if (it>=maxiter-1)
2342  {
2343  SG_SWARNING("CStatistics::fit_sigmoid(): reaching maximal iterations,"
2344  " g1=%f, g2=%f\n", g1, g2);
2345  }
2346 
2347  SG_SDEBUG("fitted sigmoid: a=%f, b=%f\n", a, b)
2348 
2349  CStatistics::SigmoidParamters result;
2350  result.a=a;
2351  result.b=b;
2352 
2353  SG_SDEBUG("leaving CStatistics::fit_sigmoid()\n")
2354  return result;
2355 }
SHOGUN Machine Learning Toolbox - Documentation