SHOGUN  v2.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 
21 #ifdef HAVE_LAPACK
22 #include <shogun/mathematics/lapack.h>
23 #endif //HAVE_LAPACK
24 using namespace shogun;
25 
26 float64_t CStatistics::mean(SGVector<float64_t> values)
27 {
28  ASSERT(values.vlen>0);
29  ASSERT(values.vector);
30 
31  float64_t sum=0;
32  for (index_t i=0; i<values.vlen; ++i)
33  sum+=values.vector[i];
34 
35  return sum/values.vlen;
36 }
37 
38 float64_t CStatistics::median(SGVector<float64_t> values, bool modify,
39  bool in_place)
40 {
41  float64_t result;
42  if (modify)
43  {
44  /* use QuickSelect method
45  * This Quickselect routine is based on the algorithm described in
46  * "Numerical recipes in C", Second Edition,
47  * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
48  * This code by Nicolas Devillard - 1998. Public domain.
49  * Adapted to SHOGUN by Heiko Strathmann
50  */
51  int32_t low;
52  int32_t high;
53  int32_t median;
54  int32_t middle;
55  int32_t l;
56  int32_t h;
57 
58  low=0;
59  high=values.vlen-1;
60  median=(low+high)/2;
61 
62  while (true)
63  {
64  if (high<=low)
65  {
66  result=values[median];
67  break;
68  }
69 
70  if (high==low+1)
71  {
72  if (values[low]>values[high])
73  CMath::CMath::swap(values[low], values[high]);
74  result=values[median];
75  break;
76  }
77 
78  middle=(low+high)/2;
79  if (values[middle]>values[high])
80  CMath::swap(values[middle], values[high]);
81  if (values[low]>values[high])
82  CMath::swap(values[low], values[high]);
83  if (values[middle]>values[low])
84  CMath::swap(values[middle], values[low]);
85 
86  CMath::swap(values[middle], values[low+1]);
87 
88  l=low+1;
89  h=high;
90  for (;;)
91  {
92  do
93  l++;
94  while (values[low]>values[l]);
95  do
96  h--;
97  while (values[h]>values[low]);
98  if (h<l)
99  break;
100  CMath::swap(values[l], values[h]);
101  }
102 
103  CMath::swap(values[low], values[h]);
104  if (h<=median)
105  low=l;
106  if (h>=median)
107  high=h-1;
108  }
109 
110  }
111  else
112  {
113  if (in_place)
114  {
115  /* use Torben method
116  * The following code is public domain.
117  * Algorithm by Torben Mogensen, implementation by N. Devillard.
118  * This code in public domain.
119  * Adapted to SHOGUN by Heiko Strathmann
120  */
121  int32_t i;
122  int32_t less;
123  int32_t greater;
124  int32_t equal;
125  float64_t min;
126  float64_t max;
127  float64_t guess;
128  float64_t maxltguess;
129  float64_t mingtguess;
130  min=max=values[0];
131  for (i=1; i<values.vlen; i++)
132  {
133  if (values[i]<min)
134  min=values[i];
135  if (values[i]>max)
136  max=values[i];
137  }
138  while (1)
139  {
140  guess=(min+max)/2;
141  less=0;
142  greater=0;
143  equal=0;
144  maxltguess=min;
145  mingtguess=max;
146  for (i=0; i<values.vlen; i++)
147  {
148  if (values[i]<guess)
149  {
150  less++;
151  if (values[i]>maxltguess)
152  maxltguess=values[i];
153  }
154  else if (values[i]>guess)
155  {
156  greater++;
157  if (values[i]<mingtguess)
158  mingtguess=values[i];
159  }
160  else
161  equal++;
162  }
163  if (less<=(values.vlen+1)/2&&greater<=(values.vlen+1)/2)
164  break;
165  else if (less>greater)
166  max=maxltguess;
167  else
168  min=mingtguess;
169  }
170 
171  if (less>=(values.vlen+1)/2)
172  result=maxltguess;
173  else if (less+equal>=(values.vlen+1)/2)
174  result=guess;
175  else
176  result=mingtguess;
177  }
178  else
179  {
180  /* copy vector and do recursive call which modifies copy */
181  SGVector<float64_t> copy(values.vlen);
182  memcpy(copy.vector, values.vector, sizeof(float64_t)*values.vlen);
183  result=median(copy, true);
184  }
185  }
186 
187  return result;
188 }
189 
190 float64_t CStatistics::matrix_median(SGMatrix<float64_t> values,
191  bool modify, bool in_place)
192 {
193  /* create a vector that uses the matrix data, dont do reference counting */
194  SGVector<float64_t> as_vector(values.matrix,
195  values.num_rows*values.num_cols, false);
196 
197  /* return vector median method */
198  return median(as_vector, modify, in_place);
199 }
200 
201 
202 float64_t CStatistics::variance(SGVector<float64_t> values)
203 {
204  ASSERT(values.vlen>1);
205  ASSERT(values.vector);
206 
207  float64_t mean=CStatistics::mean(values);
208 
209  float64_t sum_squared_diff=0;
210  for (index_t i=0; i<values.vlen; ++i)
211  sum_squared_diff+=CMath::pow(values.vector[i]-mean, 2);
212 
213  return sum_squared_diff/(values.vlen-1);
214 }
215 
216 SGVector<float64_t> CStatistics::matrix_mean(SGMatrix<float64_t> values,
217  bool col_wise)
218 {
219  ASSERT(values.num_rows>0);
220  ASSERT(values.num_cols>0);
221  ASSERT(values.matrix);
222 
223  SGVector<float64_t> result;
224 
225  if (col_wise)
226  {
227  result=SGVector<float64_t>(values.num_cols);
228  for (index_t j=0; j<values.num_cols; ++j)
229  {
230  result[j]=0;
231  for (index_t i=0; i<values.num_rows; ++i)
232  result[j]+=values(i,j);
233 
234  result[j]/=values.num_rows;
235  }
236  }
237  else
238  {
239  result=SGVector<float64_t>(values.num_rows);
240  for (index_t j=0; j<values.num_rows; ++j)
241  {
242  result[j]=0;
243  for (index_t i=0; i<values.num_cols; ++i)
244  result[j]+=values(j,i);
245 
246  result[j]/=values.num_cols;
247  }
248  }
249 
250  return result;
251 }
252 
253 SGVector<float64_t> CStatistics::matrix_variance(SGMatrix<float64_t> values,
254  bool col_wise)
255 {
256  ASSERT(values.num_rows>0);
257  ASSERT(values.num_cols>0);
258  ASSERT(values.matrix);
259 
260  /* first compute mean */
261  SGVector<float64_t> mean=CStatistics::matrix_mean(values, col_wise);
262 
263  SGVector<float64_t> result;
264 
265  if (col_wise)
266  {
267  result=SGVector<float64_t>(values.num_cols);
268  for (index_t j=0; j<values.num_cols; ++j)
269  {
270  result[j]=0;
271  for (index_t i=0; i<values.num_rows; ++i)
272  result[j]+=CMath::pow(values(i,j)-mean[j], 2);
273 
274  result[j]/=(values.num_rows-1);
275  }
276  }
277  else
278  {
279  result=SGVector<float64_t>(values.num_rows);
280  for (index_t j=0; j<values.num_rows; ++j)
281  {
282  result[j]=0;
283  for (index_t i=0; i<values.num_cols; ++i)
284  result[j]+=CMath::pow(values(j,i)-mean[j], 2);
285 
286  result[j]/=(values.num_cols-1);
287  }
288  }
289 
290  return result;
291 }
292 
293 float64_t CStatistics::std_deviation(SGVector<float64_t> values)
294 {
295  return CMath::sqrt(variance(values));
296 }
297 
298 SGVector<float64_t> CStatistics::matrix_std_deviation(
299  SGMatrix<float64_t> values, bool col_wise)
300 {
301  SGVector<float64_t> var=CStatistics::matrix_variance(values, col_wise);
302  for (index_t i=0; i<var.vlen; ++i)
303  var[i]=CMath::sqrt(var[i]);
304 
305  return var;
306 }
307 
308 #ifdef HAVE_LAPACK
309 SGMatrix<float64_t> CStatistics::covariance_matrix(
310  SGMatrix<float64_t> observations, bool in_place)
311 {
312  SGMatrix<float64_t> centered=
313  in_place ?
314  observations :
315  SGMatrix<float64_t>(observations.num_rows,
316  observations.num_cols);
317 
318  if (!in_place)
319  {
320  memcpy(centered.matrix, observations.matrix,
321  sizeof(float64_t)*observations.num_rows*observations.num_cols);
322  }
323  centered.remove_column_mean();
324 
325  /* compute 1/(m-1) * X' * X */
326  SGMatrix<float64_t> cov=SGMatrix<float64_t>::matrix_multiply(centered,
327  centered, true, false, 1.0/(observations.num_rows-1));
328 
329  return cov;
330 }
331 #endif //HAVE_LAPACK
332 
333 float64_t CStatistics::confidence_intervals_mean(SGVector<float64_t> values,
334  float64_t alpha, float64_t& conf_int_low, float64_t& conf_int_up)
335 {
336  ASSERT(values.vlen>1);
337  ASSERT(values.vector);
338 
339  /* using one sided student t distribution evaluation */
340  alpha=alpha/2;
341 
342  /* degrees of freedom */
343  int32_t deg=values.vlen-1;
344 
345  /* compute absolute value of t-value */
346  float64_t t=CMath::abs(inverse_student_t(deg, alpha));
347 
348  /* values for calculating confidence interval */
349  float64_t std_dev=CStatistics::std_deviation(values);
350  float64_t mean=CStatistics::mean(values);
351 
352  /* compute confidence interval */
353  float64_t interval=t*std_dev/CMath::sqrt((float64_t)values.vlen);
354  conf_int_low=mean-interval;
355  conf_int_up=mean+interval;
356 
357  return mean;
358 }
359 
360 float64_t CStatistics::inverse_student_t(int32_t k, float64_t p)
361 {
362  float64_t t;
363  float64_t rk;
364  float64_t z;
365  int32_t rflg;
366  float64_t result;
367 
368  if (!(k>0 && greater(p, 0)) && less(p, 1))
369  {
370  SG_SERROR("CStatistics::inverse_student_t_distribution(): "
371  "Domain error\n");
372  }
373  rk=k;
374  if (greater(p, 0.25) && less(p, 0.75))
375  {
376  if (equal(p, 0.5))
377  {
378  result=0;
379  return result;
380  }
381  z=1.0-2.0*p;
382  z=inverse_incomplete_beta(0.5, 0.5*rk, CMath::abs(z));
383  t=CMath::sqrt(rk*z/(1.0-z));
384  if (less(p, 0.5))
385  {
386  t=-t;
387  }
388  result=t;
389  return result;
390  }
391  rflg=-1;
392  if (greater_equal(p, 0.5))
393  {
394  p=1.0-p;
395  rflg=1;
396  }
397  z=inverse_incomplete_beta(0.5*rk, 0.5, 2.0*p);
398  if (less(CMath::MAX_REAL_NUMBER*z, rk))
399  {
400  result=rflg*CMath::MAX_REAL_NUMBER;
401  return result;
402  }
403  t=CMath::sqrt(rk/z-rk);
404  result=rflg*t;
405  return result;
406 }
407 
408 float64_t CStatistics::inverse_incomplete_beta(float64_t a, float64_t b,
409  float64_t y)
410 {
411  float64_t aaa;
412  float64_t bbb;
413  float64_t y0;
414  float64_t d;
415  float64_t yyy;
416  float64_t x;
417  float64_t x0;
418  float64_t x1;
419  float64_t lgm;
420  float64_t yp;
421  float64_t di;
422  float64_t dithresh;
423  float64_t yl;
424  float64_t yh;
425  float64_t xt;
426  int32_t i;
427  int32_t rflg;
428  int32_t dir;
429  int32_t nflg;
430  int32_t mainlooppos;
431  int32_t ihalve;
432  int32_t ihalvecycle;
433  int32_t newt;
434  int32_t newtcycle;
435  int32_t breaknewtcycle;
436  int32_t breakihalvecycle;
437  float64_t result;
438 
439  i=0;
440  if (!(greater_equal(y, 0) && less_equal(y, 1)))
441  {
442  SG_SERROR("CStatistics::inverse_incomplete_beta(): "
443  "Domain error\n");
444  }
445 
446  /*
447  * special cases
448  */
449  if (equal(y, 0))
450  {
451  result=0;
452  return result;
453  }
454  if (equal(y, 1.0))
455  {
456  result=1;
457  return result;
458  }
459 
460  /*
461  * these initializations are not really necessary,
462  * but without them compiler complains about 'possibly uninitialized variables'.
463  */
464  dithresh=0;
465  rflg=0;
466  aaa=0;
467  bbb=0;
468  y0=0;
469  x=0;
470  yyy=0;
471  lgm=0;
472  dir=0;
473  di=0;
474 
475  /*
476  * normal initializations
477  */
478  x0=0.0;
479  yl=0.0;
480  x1=1.0;
481  yh=1.0;
482  nflg=0;
483  mainlooppos=0;
484  ihalve=1;
485  ihalvecycle=2;
486  newt=3;
487  newtcycle=4;
488  breaknewtcycle=5;
489  breakihalvecycle=6;
490 
491  /*
492  * main loop
493  */
494  for (;;)
495  {
496 
497  /*
498  * start
499  */
500  if (mainlooppos==0)
501  {
502  if (less_equal(a, 1.0) || less_equal(b, 1.0))
503  {
504  dithresh=1.0e-6;
505  rflg=0;
506  aaa=a;
507  bbb=b;
508  y0=y;
509  x=aaa/(aaa+bbb);
510  yyy=incomplete_beta(aaa, bbb, x);
511  mainlooppos=ihalve;
512  continue;
513  }
514  else
515  {
516  dithresh=1.0e-4;
517  }
518  yp=-inverse_normal_cdf(y);
519  if (greater(y, 0.5))
520  {
521  rflg=1;
522  aaa=b;
523  bbb=a;
524  y0=1.0-y;
525  yp=-yp;
526  }
527  else
528  {
529  rflg=0;
530  aaa=a;
531  bbb=b;
532  y0=y;
533  }
534  lgm=(yp*yp-3.0)/6.0;
535  x=2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
536  d=yp*CMath::sqrt(x+lgm)/x
537  -(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))
538  *(lgm+5.0/6.0-2.0/(3.0*x));
539  d=2.0*d;
540  if (less(d, CMath::log(CMath::MIN_REAL_NUMBER)))
541  {
542  x=0;
543  break;
544  }
545  x=aaa/(aaa+bbb*CMath::exp(d));
546  yyy=incomplete_beta(aaa, bbb, x);
547  yp=(yyy-y0)/y0;
548  if (less(CMath::abs(yp), 0.2))
549  {
550  mainlooppos=newt;
551  continue;
552  }
553  mainlooppos=ihalve;
554  continue;
555  }
556 
557  /*
558  * ihalve
559  */
560  if (mainlooppos==ihalve)
561  {
562  dir=0;
563  di=0.5;
564  i=0;
565  mainlooppos=ihalvecycle;
566  continue;
567  }
568 
569  /*
570  * ihalvecycle
571  */
572  if (mainlooppos==ihalvecycle)
573  {
574  if (i<=99)
575  {
576  if (i!=0)
577  {
578  x=x0+di*(x1-x0);
579  if (equal(x, 1.0))
580  {
581  x=1.0-CMath::MACHINE_EPSILON;
582  }
583  if (equal(x, 0.0))
584  {
585  di=0.5;
586  x=x0+di*(x1-x0);
587  if (equal(x, 0.0))
588  {
589  break;
590  }
591  }
592  yyy=incomplete_beta(aaa, bbb, x);
593  yp=(x1-x0)/(x1+x0);
594  if (less(CMath::abs(yp), dithresh))
595  {
596  mainlooppos=newt;
597  continue;
598  }
599  yp=(yyy-y0)/y0;
600  if (less(CMath::abs(yp), dithresh))
601  {
602  mainlooppos=newt;
603  continue;
604  }
605  }
606  if (less(yyy, y0))
607  {
608  x0=x;
609  yl=yyy;
610  if (dir<0)
611  {
612  dir=0;
613  di=0.5;
614  }
615  else
616  {
617  if (dir>3)
618  {
619  di=1.0-(1.0-di)*(1.0-di);
620  }
621  else
622  {
623  if (dir>1)
624  {
625  di=0.5*di+0.5;
626  }
627  else
628  {
629  di=(y0-yyy)/(yh-yl);
630  }
631  }
632  }
633  dir=dir+1;
634  if (greater(x0, 0.75))
635  {
636  if (rflg==1)
637  {
638  rflg=0;
639  aaa=a;
640  bbb=b;
641  y0=y;
642  }
643  else
644  {
645  rflg=1;
646  aaa=b;
647  bbb=a;
648  y0=1.0-y;
649  }
650  x=1.0-x;
651  yyy=incomplete_beta(aaa, bbb, x);
652  x0=0.0;
653  yl=0.0;
654  x1=1.0;
655  yh=1.0;
656  mainlooppos=ihalve;
657  continue;
658  }
659  }
660  else
661  {
662  x1=x;
663  if (rflg==1 && less(x1, CMath::MACHINE_EPSILON))
664  {
665  x=0.0;
666  break;
667  }
668  yh=yyy;
669  if (dir>0)
670  {
671  dir=0;
672  di=0.5;
673  }
674  else
675  {
676  if (dir<-3)
677  {
678  di=di*di;
679  }
680  else
681  {
682  if (dir<-1)
683  {
684  di=0.5*di;
685  }
686  else
687  {
688  di=(yyy-y0)/(yh-yl);
689  }
690  }
691  }
692  dir=dir-1;
693  }
694  i=i+1;
695  mainlooppos=ihalvecycle;
696  continue;
697  }
698  else
699  {
700  mainlooppos=breakihalvecycle;
701  continue;
702  }
703  }
704 
705  /*
706  * breakihalvecycle
707  */
708  if (mainlooppos==breakihalvecycle)
709  {
710  if (greater_equal(x0, 1.0))
711  {
712  x=1.0-CMath::MACHINE_EPSILON;
713  break;
714  }
715  if (less_equal(x, 0.0))
716  {
717  x=0.0;
718  break;
719  }
720  mainlooppos=newt;
721  continue;
722  }
723 
724  /*
725  * newt
726  */
727  if (mainlooppos==newt)
728  {
729  if (nflg!=0)
730  {
731  break;
732  }
733  nflg=1;
734  lgm=lgamma(aaa+bbb)-lgamma(aaa)-lgamma(bbb);
735  i=0;
736  mainlooppos=newtcycle;
737  continue;
738  }
739 
740  /*
741  * newtcycle
742  */
743  if (mainlooppos==newtcycle)
744  {
745  if (i<=7)
746  {
747  if (i!=0)
748  {
749  yyy=incomplete_beta(aaa, bbb, x);
750  }
751  if (less(yyy, yl))
752  {
753  x=x0;
754  yyy=yl;
755  }
756  else
757  {
758  if (greater(yyy, yh))
759  {
760  x=x1;
761  yyy=yh;
762  }
763  else
764  {
765  if (less(yyy, y0))
766  {
767  x0=x;
768  yl=yyy;
769  }
770  else
771  {
772  x1=x;
773  yh=yyy;
774  }
775  }
776  }
777  if (equal(x, 1.0) || equal(x, 0.0))
778  {
779  mainlooppos=breaknewtcycle;
780  continue;
781  }
782  d=(aaa-1.0)*CMath::log(x)+(bbb-1.0)*CMath::log(1.0-x)+lgm;
783  if (less(d, CMath::log(CMath::MIN_REAL_NUMBER)))
784  {
785  break;
786  }
787  if (greater(d, CMath::log(CMath::MAX_REAL_NUMBER)))
788  {
789  mainlooppos=breaknewtcycle;
790  continue;
791  }
792  d=CMath::exp(d);
793  d=(yyy-y0)/d;
794  xt=x-d;
795  if (less_equal(xt, x0))
796  {
797  yyy=(x-x0)/(x1-x0);
798  xt=x0+0.5*yyy*(x-x0);
799  if (less_equal(xt, 0.0))
800  {
801  mainlooppos=breaknewtcycle;
802  continue;
803  }
804  }
805  if (greater_equal(xt, x1))
806  {
807  yyy=(x1-x)/(x1-x0);
808  xt=x1-0.5*yyy*(x1-x);
809  if (greater_equal(xt, 1.0))
810  {
811  mainlooppos=breaknewtcycle;
812  continue;
813  }
814  }
815  x=xt;
816  if (less(CMath::abs(d/x), 128.0*CMath::MACHINE_EPSILON))
817  {
818  break;
819  }
820  i=i+1;
821  mainlooppos=newtcycle;
822  continue;
823  }
824  else
825  {
826  mainlooppos=breaknewtcycle;
827  continue;
828  }
829  }
830 
831  /*
832  * breaknewtcycle
833  */
834  if (mainlooppos==breaknewtcycle)
835  {
836  dithresh=256.0*CMath::MACHINE_EPSILON;
837  mainlooppos=ihalve;
838  continue;
839  }
840  }
841 
842  /*
843  * done
844  */
845  if (rflg!=0)
846  {
847  if (less_equal(x, CMath::MACHINE_EPSILON))
848  {
849  x=1.0-CMath::MACHINE_EPSILON;
850  }
851  else
852  {
853  x=1.0-x;
854  }
855  }
856  result=x;
857  return result;
858 }
859 
860 float64_t CStatistics::incomplete_beta(float64_t a, float64_t b, float64_t x)
861 {
862  float64_t t;
863  float64_t xc;
864  float64_t w;
865  float64_t y;
866  int32_t flag;
867  float64_t big;
868  float64_t biginv;
869  float64_t maxgam;
870  float64_t minlog;
871  float64_t maxlog;
872  float64_t result;
873 
874  big=4.503599627370496e15;
875  biginv=2.22044604925031308085e-16;
876  maxgam=171.624376956302725;
877  minlog=CMath::log(CMath::MIN_REAL_NUMBER);
878  maxlog=CMath::log(CMath::MAX_REAL_NUMBER);
879 
880  if (!(greater(a, 0) && greater(b, 0)))
881  {
882  SG_SERROR("CStatistics::incomplete_beta(): "
883  "Domain error\n");
884  }
885 
886  if (!(greater_equal(x, 0) && less_equal(x, 1)))
887  {
888  SG_SERROR("CStatistics::incomplete_beta(): "
889  "Domain error\n");
890  }
891 
892  if (equal(x, 0))
893  {
894  result=0;
895  return result;
896  }
897  if (equal(x, 1))
898  {
899  result=1;
900  return result;
901  }
902  flag=0;
903  if (less_equal(b*x, 1.0) && less_equal(x, 0.95))
904  {
905  result=ibetaf_incompletebetaps(a, b, x, maxgam);
906  return result;
907  }
908  w=1.0-x;
909  if (greater(x, a/(a+b)))
910  {
911  flag=1;
912  t=a;
913  a=b;
914  b=t;
915  xc=x;
916  x=w;
917  }
918  else
919  {
920  xc=w;
921  }
922  if ((flag==1 && less_equal(b*x, 1.0)) && less_equal(x, 0.95))
923  {
924  t=ibetaf_incompletebetaps(a, b, x, maxgam);
925  if (less_equal(t, CMath::MACHINE_EPSILON))
926  {
927  result=1.0-CMath::MACHINE_EPSILON;
928  }
929  else
930  {
931  result=1.0-t;
932  }
933  return result;
934  }
935  y=x*(a+b-2.0)-(a-1.0);
936  if (less(y, 0.0))
937  {
938  w=ibetaf_incompletebetafe(a, b, x, big, biginv);
939  }
940  else
941  {
942  w=ibetaf_incompletebetafe2(a, b, x, big, biginv)/xc;
943  }
944  y=a*CMath::log(x);
945  t=b*CMath::log(xc);
946  if ((less(a+b, maxgam) && less(CMath::abs(y), maxlog))
947  && less(CMath::abs(t), maxlog))
948  {
949  t=CMath::pow(xc, b);
950  t=t*CMath::pow(x, a);
951  t=t/a;
952  t=t*w;
953  t=t*(tgamma(a+b)/(tgamma(a)*tgamma(b)));
954  if (flag==1)
955  {
956  if (less_equal(t, CMath::MACHINE_EPSILON))
957  {
958  result=1.0-CMath::MACHINE_EPSILON;
959  }
960  else
961  {
962  result=1.0-t;
963  }
964  }
965  else
966  {
967  result=t;
968  }
969  return result;
970  }
971  y=y+t+lgamma(a+b)-lgamma(a)-lgamma(b);
972  y=y+CMath::log(w/a);
973  if (less(y, minlog))
974  {
975  t=0.0;
976  }
977  else
978  {
979  t=CMath::exp(y);
980  }
981  if (flag==1)
982  {
983  if (less_equal(t, CMath::MACHINE_EPSILON))
984  {
985  t=1.0-CMath::MACHINE_EPSILON;
986  }
987  else
988  {
989  t=1.0-t;
990  }
991  }
992  result=t;
993  return result;
994 }
995 
996 float64_t CStatistics::inverse_normal_cdf(float64_t y0, float64_t mean,
997  float64_t std_dev)
998 {
999  return inverse_normal_cdf(y0)*std_dev+mean;
1000 }
1001 
1002 float64_t CStatistics::inverse_normal_cdf(float64_t y0)
1003 {
1004  float64_t expm2;
1005  float64_t s2pi;
1006  float64_t x;
1007  float64_t y;
1008  float64_t z;
1009  float64_t y2;
1010  float64_t x0;
1011  float64_t x1;
1012  int32_t code;
1013  float64_t p0;
1014  float64_t q0;
1015  float64_t p1;
1016  float64_t q1;
1017  float64_t p2;
1018  float64_t q2;
1019  float64_t result;
1020 
1021  expm2=0.13533528323661269189;
1022  s2pi=2.50662827463100050242;
1023  if (less_equal(y0, 0))
1024  {
1025  result=-CMath::MAX_REAL_NUMBER;
1026  return result;
1027  }
1028  if (greater_equal(y0, 1))
1029  {
1030  result=CMath::MAX_REAL_NUMBER;
1031  return result;
1032  }
1033  code=1;
1034  y=y0;
1035  if (greater(y, 1.0-expm2))
1036  {
1037  y=1.0-y;
1038  code=0;
1039  }
1040  if (greater(y, expm2))
1041  {
1042  y=y-0.5;
1043  y2=y*y;
1044  p0=-59.9633501014107895267;
1045  p0=98.0010754185999661536+y2*p0;
1046  p0=-56.6762857469070293439+y2*p0;
1047  p0=13.9312609387279679503+y2*p0;
1048  p0=-1.23916583867381258016+y2*p0;
1049  q0=1;
1050  q0=1.95448858338141759834+y2*q0;
1051  q0=4.67627912898881538453+y2*q0;
1052  q0=86.3602421390890590575+y2*q0;
1053  q0=-225.462687854119370527+y2*q0;
1054  q0=200.260212380060660359+y2*q0;
1055  q0=-82.0372256168333339912+y2*q0;
1056  q0=15.9056225126211695515+y2*q0;
1057  q0=-1.18331621121330003142+y2*q0;
1058  x=y+y*y2*p0/q0;
1059  x=x*s2pi;
1060  result=x;
1061  return result;
1062  }
1063  x=CMath::sqrt(-2.0*CMath::log(y));
1064  x0=x-CMath::log(x)/x;
1065  z=1.0/x;
1066  if (less(x, 8.0))
1067  {
1068  p1=4.05544892305962419923;
1069  p1=31.5251094599893866154+z*p1;
1070  p1=57.1628192246421288162+z*p1;
1071  p1=44.0805073893200834700+z*p1;
1072  p1=14.6849561928858024014+z*p1;
1073  p1=2.18663306850790267539+z*p1;
1074  p1=-1.40256079171354495875*0.1+z*p1;
1075  p1=-3.50424626827848203418*0.01+z*p1;
1076  p1=-8.57456785154685413611*0.0001+z*p1;
1077  q1=1;
1078  q1=15.7799883256466749731+z*q1;
1079  q1=45.3907635128879210584+z*q1;
1080  q1=41.3172038254672030440+z*q1;
1081  q1=15.0425385692907503408+z*q1;
1082  q1=2.50464946208309415979+z*q1;
1083  q1=-1.42182922854787788574*0.1+z*q1;
1084  q1=-3.80806407691578277194*0.01+z*q1;
1085  q1=-9.33259480895457427372*0.0001+z*q1;
1086  x1=z*p1/q1;
1087  }
1088  else
1089  {
1090  p2=3.23774891776946035970;
1091  p2=6.91522889068984211695+z*p2;
1092  p2=3.93881025292474443415+z*p2;
1093  p2=1.33303460815807542389+z*p2;
1094  p2=2.01485389549179081538*0.1+z*p2;
1095  p2=1.23716634817820021358*0.01+z*p2;
1096  p2=3.01581553508235416007*0.0001+z*p2;
1097  p2=2.65806974686737550832*0.000001+z*p2;
1098  p2=6.23974539184983293730*0.000000001+z*p2;
1099  q2=1;
1100  q2=6.02427039364742014255+z*q2;
1101  q2=3.67983563856160859403+z*q2;
1102  q2=1.37702099489081330271+z*q2;
1103  q2=2.16236993594496635890*0.1+z*q2;
1104  q2=1.34204006088543189037*0.01+z*q2;
1105  q2=3.28014464682127739104*0.0001+z*q2;
1106  q2=2.89247864745380683936*0.000001+z*q2;
1107  q2=6.79019408009981274425*0.000000001+z*q2;
1108  x1=z*p2/q2;
1109  }
1110  x=x0-x1;
1111  if (code!=0)
1112  {
1113  x=-x;
1114  }
1115  result=x;
1116  return result;
1117 }
1118 
1119 float64_t CStatistics::ibetaf_incompletebetaps(float64_t a, float64_t b,
1120  float64_t x, float64_t maxgam)
1121 {
1122  float64_t s;
1123  float64_t t;
1124  float64_t u;
1125  float64_t v;
1126  float64_t n;
1127  float64_t t1;
1128  float64_t z;
1129  float64_t ai;
1130  float64_t result;
1131 
1132  ai=1.0/a;
1133  u=(1.0-b)*x;
1134  v=u/(a+1.0);
1135  t1=v;
1136  t=u;
1137  n=2.0;
1138  s=0.0;
1139  z=CMath::MACHINE_EPSILON*ai;
1140  while (greater(CMath::abs(v), z))
1141  {
1142  u=(n-b)*x/n;
1143  t=t*u;
1144  v=t/(a+n);
1145  s=s+v;
1146  n=n+1.0;
1147  }
1148  s=s+t1;
1149  s=s+ai;
1150  u=a*CMath::log(x);
1151  if (less(a+b, maxgam)
1152  && less(CMath::abs(u), CMath::log(CMath::MAX_REAL_NUMBER)))
1153  {
1154  t=tgamma(a+b)/(tgamma(a)*tgamma(b));
1155  s=s*t*CMath::pow(x, a);
1156  }
1157  else
1158  {
1159  t=lgamma(a+b)-lgamma(a)-lgamma(b)+u+CMath::log(s);
1160  if (less(t, CMath::log(CMath::MIN_REAL_NUMBER)))
1161  {
1162  s=0.0;
1163  }
1164  else
1165  {
1166  s=CMath::exp(t);
1167  }
1168  }
1169  result=s;
1170  return result;
1171 }
1172 
1173 float64_t CStatistics::ibetaf_incompletebetafe(float64_t a, float64_t b,
1174  float64_t x, float64_t big, float64_t biginv)
1175 {
1176  float64_t xk;
1177  float64_t pk;
1178  float64_t pkm1;
1179  float64_t pkm2;
1180  float64_t qk;
1181  float64_t qkm1;
1182  float64_t qkm2;
1183  float64_t k1;
1184  float64_t k2;
1185  float64_t k3;
1186  float64_t k4;
1187  float64_t k5;
1188  float64_t k6;
1189  float64_t k7;
1190  float64_t k8;
1191  float64_t r;
1192  float64_t t;
1193  float64_t ans;
1194  float64_t thresh;
1195  int32_t n;
1196  float64_t result;
1197 
1198  k1=a;
1199  k2=a+b;
1200  k3=a;
1201  k4=a+1.0;
1202  k5=1.0;
1203  k6=b-1.0;
1204  k7=k4;
1205  k8=a+2.0;
1206  pkm2=0.0;
1207  qkm2=1.0;
1208  pkm1=1.0;
1209  qkm1=1.0;
1210  ans=1.0;
1211  r=1.0;
1212  n=0;
1213  thresh=3.0*CMath::MACHINE_EPSILON;
1214  do
1215  {
1216  xk=-x*k1*k2/(k3*k4);
1217  pk=pkm1+pkm2*xk;
1218  qk=qkm1+qkm2*xk;
1219  pkm2=pkm1;
1220  pkm1=pk;
1221  qkm2=qkm1;
1222  qkm1=qk;
1223  xk=x*k5*k6/(k7*k8);
1224  pk=pkm1+pkm2*xk;
1225  qk=qkm1+qkm2*xk;
1226  pkm2=pkm1;
1227  pkm1=pk;
1228  qkm2=qkm1;
1229  qkm1=qk;
1230  if (not_equal(qk, 0))
1231  {
1232  r=pk/qk;
1233  }
1234  if (not_equal(r, 0))
1235  {
1236  t=CMath::abs((ans-r)/r);
1237  ans=r;
1238  }
1239  else
1240  {
1241  t=1.0;
1242  }
1243  if (less(t, thresh))
1244  {
1245  break;
1246  }
1247  k1=k1+1.0;
1248  k2=k2+1.0;
1249  k3=k3+2.0;
1250  k4=k4+2.0;
1251  k5=k5+1.0;
1252  k6=k6-1.0;
1253  k7=k7+2.0;
1254  k8=k8+2.0;
1255  if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1256  {
1257  pkm2=pkm2*biginv;
1258  pkm1=pkm1*biginv;
1259  qkm2=qkm2*biginv;
1260  qkm1=qkm1*biginv;
1261  }
1262  if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1263  {
1264  pkm2=pkm2*big;
1265  pkm1=pkm1*big;
1266  qkm2=qkm2*big;
1267  qkm1=qkm1*big;
1268  }
1269  n=n+1;
1270  }
1271  while (n!=300);
1272  result=ans;
1273  return result;
1274 }
1275 
1276 float64_t CStatistics::ibetaf_incompletebetafe2(float64_t a, float64_t b,
1277  float64_t x, float64_t big, float64_t biginv)
1278 {
1279  float64_t xk;
1280  float64_t pk;
1281  float64_t pkm1;
1282  float64_t pkm2;
1283  float64_t qk;
1284  float64_t qkm1;
1285  float64_t qkm2;
1286  float64_t k1;
1287  float64_t k2;
1288  float64_t k3;
1289  float64_t k4;
1290  float64_t k5;
1291  float64_t k6;
1292  float64_t k7;
1293  float64_t k8;
1294  float64_t r;
1295  float64_t t;
1296  float64_t ans;
1297  float64_t z;
1298  float64_t thresh;
1299  int32_t n;
1300  float64_t result;
1301 
1302  k1=a;
1303  k2=b-1.0;
1304  k3=a;
1305  k4=a+1.0;
1306  k5=1.0;
1307  k6=a+b;
1308  k7=a+1.0;
1309  k8=a+2.0;
1310  pkm2=0.0;
1311  qkm2=1.0;
1312  pkm1=1.0;
1313  qkm1=1.0;
1314  z=x/(1.0-x);
1315  ans=1.0;
1316  r=1.0;
1317  n=0;
1318  thresh=3.0*CMath::MACHINE_EPSILON;
1319  do
1320  {
1321  xk=-z*k1*k2/(k3*k4);
1322  pk=pkm1+pkm2*xk;
1323  qk=qkm1+qkm2*xk;
1324  pkm2=pkm1;
1325  pkm1=pk;
1326  qkm2=qkm1;
1327  qkm1=qk;
1328  xk=z*k5*k6/(k7*k8);
1329  pk=pkm1+pkm2*xk;
1330  qk=qkm1+qkm2*xk;
1331  pkm2=pkm1;
1332  pkm1=pk;
1333  qkm2=qkm1;
1334  qkm1=qk;
1335  if (not_equal(qk, 0))
1336  {
1337  r=pk/qk;
1338  }
1339  if (not_equal(r, 0))
1340  {
1341  t=CMath::abs((ans-r)/r);
1342  ans=r;
1343  }
1344  else
1345  {
1346  t=1.0;
1347  }
1348  if (less(t, thresh))
1349  {
1350  break;
1351  }
1352  k1=k1+1.0;
1353  k2=k2-1.0;
1354  k3=k3+2.0;
1355  k4=k4+2.0;
1356  k5=k5+1.0;
1357  k6=k6+1.0;
1358  k7=k7+2.0;
1359  k8=k8+2.0;
1360  if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1361  {
1362  pkm2=pkm2*biginv;
1363  pkm1=pkm1*biginv;
1364  qkm2=qkm2*biginv;
1365  qkm1=qkm1*biginv;
1366  }
1367  if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1368  {
1369  pkm2=pkm2*big;
1370  pkm1=pkm1*big;
1371  qkm2=qkm2*big;
1372  qkm1=qkm1*big;
1373  }
1374  n=n+1;
1375  }
1376  while (n!=300);
1377  result=ans;
1378  return result;
1379 }
1380 
1381 float64_t CStatistics::incomplete_gamma(float64_t a, float64_t x)
1382 {
1383  float64_t igammaepsilon;
1384  float64_t ans;
1385  float64_t ax;
1386  float64_t c;
1387  float64_t r;
1388  float64_t result;
1389 
1390  igammaepsilon=0.000000000000001;
1391  if (less_equal(x, 0) || less_equal(a, 0))
1392  {
1393  result=0;
1394  return result;
1395  }
1396  if (greater(x, 1) && greater(x, a))
1397  {
1398  result=1-incomplete_gamma_completed(a, x);
1399  return result;
1400  }
1401  ax=a*CMath::log(x)-x-lgamma(a);
1402  if (less(ax, -709.78271289338399))
1403  {
1404  result=0;
1405  return result;
1406  }
1407  ax=CMath::exp(ax);
1408  r=a;
1409  c=1;
1410  ans=1;
1411  do
1412  {
1413  r=r+1;
1414  c=c*x/r;
1415  ans=ans+c;
1416  }
1417  while (greater(c/ans, igammaepsilon));
1418  result=ans*ax/a;
1419  return result;
1420 }
1421 
1422 float64_t CStatistics::incomplete_gamma_completed(float64_t a, float64_t x)
1423 {
1424  float64_t igammaepsilon;
1425  float64_t igammabignumber;
1426  float64_t igammabignumberinv;
1427  float64_t ans;
1428  float64_t ax;
1429  float64_t c;
1430  float64_t yc;
1431  float64_t r;
1432  float64_t t;
1433  float64_t y;
1434  float64_t z;
1435  float64_t pk;
1436  float64_t pkm1;
1437  float64_t pkm2;
1438  float64_t qk;
1439  float64_t qkm1;
1440  float64_t qkm2;
1441  float64_t result;
1442 
1443  igammaepsilon=0.000000000000001;
1444  igammabignumber=4503599627370496.0;
1445  igammabignumberinv=2.22044604925031308085*0.0000000000000001;
1446  if (less_equal(x, 0) || less_equal(a, 0))
1447  {
1448  result=1;
1449  return result;
1450  }
1451  if (less(x, 1) || less(x, a))
1452  {
1453  result=1-incomplete_gamma(a, x);
1454  return result;
1455  }
1456  ax=a*CMath::log(x)-x-lgamma(a);
1457  if (less(ax, -709.78271289338399))
1458  {
1459  result=0;
1460  return result;
1461  }
1462  ax=CMath::exp(ax);
1463  y=1-a;
1464  z=x+y+1;
1465  c=0;
1466  pkm2=1;
1467  qkm2=x;
1468  pkm1=x+1;
1469  qkm1=z*x;
1470  ans=pkm1/qkm1;
1471  do
1472  {
1473  c=c+1;
1474  y=y+1;
1475  z=z+2;
1476  yc=y*c;
1477  pk=pkm1*z-pkm2*yc;
1478  qk=qkm1*z-qkm2*yc;
1479  if (not_equal(qk, 0))
1480  {
1481  r=pk/qk;
1482  t=CMath::abs((ans-r)/r);
1483  ans=r;
1484  }
1485  else
1486  {
1487  t=1;
1488  }
1489  pkm2=pkm1;
1490  pkm1=pk;
1491  qkm2=qkm1;
1492  qkm1=qk;
1493  if (greater(CMath::abs(pk), igammabignumber))
1494  {
1495  pkm2=pkm2*igammabignumberinv;
1496  pkm1=pkm1*igammabignumberinv;
1497  qkm2=qkm2*igammabignumberinv;
1498  qkm1=qkm1*igammabignumberinv;
1499  }
1500  }
1501  while (greater(t, igammaepsilon));
1502  result=ans*ax;
1503  return result;
1504 }
1505 
1506 float64_t CStatistics::gamma_cdf(float64_t x, float64_t a, float64_t b)
1507 {
1508  /* definition of wikipedia: incomplete gamma devised by true gamma */
1509  return incomplete_gamma(a, x/b);
1510 }
1511 
1512 float64_t CStatistics::inverse_gamma_cdf(float64_t p, float64_t a,
1513  float64_t b)
1514 {
1515  /* inverse of gamma(a,b) CDF is
1516  * inverse_incomplete_gamma_completed(a, 1. - p) * b */
1517  return inverse_incomplete_gamma_completed(a, 1-p)*b;
1518 }
1519 
1520 float64_t CStatistics::inverse_incomplete_gamma_completed(float64_t a,
1521  float64_t y0)
1522 {
1523  float64_t igammaepsilon;
1524  float64_t iinvgammabignumber;
1525  float64_t x0;
1526  float64_t x1;
1527  float64_t x;
1528  float64_t yl;
1529  float64_t yh;
1530  float64_t y;
1531  float64_t d;
1532  float64_t lgm;
1533  float64_t dithresh;
1534  int32_t i;
1535  int32_t dir;
1536  float64_t result;
1537 
1538  igammaepsilon=0.000000000000001;
1539  iinvgammabignumber=4503599627370496.0;
1540  x0=iinvgammabignumber;
1541  yl=0;
1542  x1=0;
1543  yh=1;
1544  dithresh=5*igammaepsilon;
1545  d=1/(9*a);
1546  y=1-d-inverse_normal_cdf(y0)*CMath::sqrt(d);
1547  x=a*y*y*y;
1548  lgm=lgamma(a);
1549  i=0;
1550  while (i<10)
1551  {
1552  if (greater(x, x0) || less(x, x1))
1553  {
1554  d=0.0625;
1555  break;
1556  }
1557  y=incomplete_gamma_completed(a, x);
1558  if (less(y, yl) || greater(y, yh))
1559  {
1560  d=0.0625;
1561  break;
1562  }
1563  if (less(y, y0))
1564  {
1565  x0=x;
1566  yl=y;
1567  }
1568  else
1569  {
1570  x1=x;
1571  yh=y;
1572  }
1573  d=(a-1)*CMath::log(x)-x-lgm;
1574  if (less(d, -709.78271289338399))
1575  {
1576  d=0.0625;
1577  break;
1578  }
1579  d=-CMath::exp(d);
1580  d=(y-y0)/d;
1581  if (less(CMath::abs(d/x), igammaepsilon))
1582  {
1583  result=x;
1584  return result;
1585  }
1586  x=x-d;
1587  i=i+1;
1588  }
1589  if (equal(x0, iinvgammabignumber))
1590  {
1591  if (less_equal(x, 0))
1592  {
1593  x=1;
1594  }
1595  while (equal(x0, iinvgammabignumber))
1596  {
1597  x=(1+d)*x;
1598  y=incomplete_gamma_completed(a, x);
1599  if (less(y, y0))
1600  {
1601  x0=x;
1602  yl=y;
1603  break;
1604  }
1605  d=d+d;
1606  }
1607  }
1608  d=0.5;
1609  dir=0;
1610  i=0;
1611  while (i<400)
1612  {
1613  x=x1+d*(x0-x1);
1614  y=incomplete_gamma_completed(a, x);
1615  lgm=(x0-x1)/(x1+x0);
1616  if (less(CMath::abs(lgm), dithresh))
1617  {
1618  break;
1619  }
1620  lgm=(y-y0)/y0;
1621  if (less(CMath::abs(lgm), dithresh))
1622  {
1623  break;
1624  }
1625  if (less_equal(x, 0.0))
1626  {
1627  break;
1628  }
1629  if (greater_equal(y, y0))
1630  {
1631  x1=x;
1632  yh=y;
1633  if (dir<0)
1634  {
1635  dir=0;
1636  d=0.5;
1637  }
1638  else
1639  {
1640  if (dir>1)
1641  {
1642  d=0.5*d+0.5;
1643  }
1644  else
1645  {
1646  d=(y0-yl)/(yh-yl);
1647  }
1648  }
1649  dir=dir+1;
1650  }
1651  else
1652  {
1653  x0=x;
1654  yl=y;
1655  if (dir>0)
1656  {
1657  dir=0;
1658  d=0.5;
1659  }
1660  else
1661  {
1662  if (dir<-1)
1663  {
1664  d=0.5*d;
1665  }
1666  else
1667  {
1668  d=(y0-yl)/(yh-yl);
1669  }
1670  }
1671  dir=dir-1;
1672  }
1673  i=i+1;
1674  }
1675  result=x;
1676  return result;
1677 }
1678 
1679 float64_t CStatistics::normal_cdf(float64_t x, float64_t std_dev)
1680 {
1681  return 0.5*(error_function(x/std_dev/1.41421356237309504880)+1);
1682 }
1683 
1684 float64_t CStatistics::error_function(float64_t x)
1685 {
1686  float64_t xsq;
1687  float64_t s;
1688  float64_t p;
1689  float64_t q;
1690  float64_t result;
1691 
1692  s=CMath::sign(x);
1693  x=CMath::abs(x);
1694  if (less(x, 0.5))
1695  {
1696  xsq=x*x;
1697  p=0.007547728033418631287834;
1698  p=0.288805137207594084924010+xsq*p;
1699  p=14.3383842191748205576712+xsq*p;
1700  p=38.0140318123903008244444+xsq*p;
1701  p=3017.82788536507577809226+xsq*p;
1702  p=7404.07142710151470082064+xsq*p;
1703  p=80437.3630960840172832162+xsq*p;
1704  q=0.0;
1705  q=1.00000000000000000000000+xsq*q;
1706  q=38.0190713951939403753468+xsq*q;
1707  q=658.070155459240506326937+xsq*q;
1708  q=6379.60017324428279487120+xsq*q;
1709  q=34216.5257924628539769006+xsq*q;
1710  q=80437.3630960840172826266+xsq*q;
1711  result=s*1.1283791670955125738961589031*x*p/q;
1712  return result;
1713  }
1714  if (greater_equal(x, 10))
1715  {
1716  result=s;
1717  return result;
1718  }
1719  result=s*(1-error_function_complement(x));
1720  return result;
1721 }
1722 
1723 float64_t CStatistics::error_function_complement(float64_t x)
1724 {
1725  float64_t p;
1726  float64_t q;
1727  float64_t result;
1728 
1729  if (less(x, 0))
1730  {
1731  result=2-error_function_complement(-x);
1732  return result;
1733  }
1734  if (less(x, 0.5))
1735  {
1736  result=1.0-error_function(x);
1737  return result;
1738  }
1739  if (greater_equal(x, 10))
1740  {
1741  result=0;
1742  return result;
1743  }
1744  p=0.0;
1745  p=0.5641877825507397413087057563+x*p;
1746  p=9.675807882987265400604202961+x*p;
1747  p=77.08161730368428609781633646+x*p;
1748  p=368.5196154710010637133875746+x*p;
1749  p=1143.262070703886173606073338+x*p;
1750  p=2320.439590251635247384768711+x*p;
1751  p=2898.0293292167655611275846+x*p;
1752  p=1826.3348842295112592168999+x*p;
1753  q=1.0;
1754  q=17.14980943627607849376131193+x*q;
1755  q=137.1255960500622202878443578+x*q;
1756  q=661.7361207107653469211984771+x*q;
1757  q=2094.384367789539593790281779+x*q;
1758  q=4429.612803883682726711528526+x*q;
1759  q=6089.5424232724435504633068+x*q;
1760  q=4958.82756472114071495438422+x*q;
1761  q=1826.3348842295112595576438+x*q;
1762  result=CMath::exp(-x*x)*p/q;
1763  return result;
1764 }
1765 
1766 SGVector<float64_t> CStatistics::fishers_exact_test_for_multiple_2x3_tables(
1767  SGMatrix<float64_t> tables)
1768 {
1769  SGMatrix<float64_t> table(NULL, 2, 3, false);
1770  int32_t len=tables.num_cols/3;
1771 
1772  SGVector<float64_t> v(len);
1773  for (int32_t i=0; i<len; i++)
1774  {
1775  table.matrix=&tables.matrix[2*3*i];
1776  v.vector[i]=fishers_exact_test_for_2x3_table(table);
1777  }
1778  return v;
1779 }
1780 
1781 float64_t CStatistics::fishers_exact_test_for_2x3_table(
1782  SGMatrix<float64_t> table)
1783 {
1784  ASSERT(table.num_rows==2);
1785  ASSERT(table.num_cols==3);
1786 
1787  int32_t m_len=3+2;
1788  float64_t* m=SG_MALLOC(float64_t, 3+2);
1789  m[0]=table.matrix[0]+table.matrix[2]+table.matrix[4];
1790  m[1]=table.matrix[1]+table.matrix[3]+table.matrix[5];
1791  m[2]=table.matrix[0]+table.matrix[1];
1792  m[3]=table.matrix[2]+table.matrix[3];
1793  m[4]=table.matrix[4]+table.matrix[5];
1794 
1795  float64_t n=SGVector<float64_t>::sum(m, m_len)/2.0;
1796  int32_t x_len=2*3*CMath::sq(SGVector<float64_t>::max(m, m_len));
1797  float64_t* x=SG_MALLOC(float64_t, x_len);
1798  SGVector<float64_t>::fill_vector(x, x_len, 0.0);
1799 
1800  float64_t log_nom=0.0;
1801  for (int32_t i=0; i<3+2; i++)
1802  log_nom+=lgamma(m[i]+1);
1803  log_nom-=lgamma(n+1.0);
1804 
1805  float64_t log_denomf=0;
1806  floatmax_t log_denom=0;
1807 
1808  for (int32_t i=0; i<3*2; i++)
1809  {
1810  log_denom+=lgammal((floatmax_t)table.matrix[i]+1);
1811  log_denomf+=lgammal((floatmax_t)table.matrix[i]+1);
1812  }
1813 
1814  floatmax_t prob_table_log=log_nom-log_denom;
1815 
1816  int32_t dim1=CMath::min(m[0], m[2]);
1817 
1818  //traverse all possible tables with given m
1819  int32_t counter=0;
1820  for (int32_t k=0; k<=dim1; k++)
1821  {
1822  for (int32_t l=CMath::max(0.0, m[0]-m[4]-k);
1823  l<=CMath::min(m[0]-k, m[3]); l++)
1824  {
1825  x[0+0*2+counter*2*3]=k;
1826  x[0+1*2+counter*2*3]=l;
1827  x[0+2*2+counter*2*3]=m[0]-x[0+0*2+counter*2*3]-x[0+1*2+counter*2*3];
1828  x[1+0*2+counter*2*3]=m[2]-x[0+0*2+counter*2*3];
1829  x[1+1*2+counter*2*3]=m[3]-x[0+1*2+counter*2*3];
1830  x[1+2*2+counter*2*3]=m[4]-x[0+2*2+counter*2*3];
1831 
1832  counter++;
1833  }
1834  }
1835 
1836 //#define DEBUG_FISHER_TABLE
1837 #ifdef DEBUG_FISHER_TABLE
1838  SG_SPRINT("counter=%d\n", counter);
1839  SG_SPRINT("dim1=%d\n", dim1);
1840  SG_SPRINT("l=%g...%g\n", CMath::max(0.0,m[0]-m[4]-0), CMath::min(m[0]-0, m[3]));
1841  SG_SPRINT("n=%g\n", n);
1842  SG_SPRINT("prob_table_log=%.18Lg\n", prob_table_log);
1843  SG_SPRINT("log_denomf=%.18g\n", log_denomf);
1844  SG_SPRINT("log_denom=%.18Lg\n", log_denom);
1845  SG_SPRINT("log_nom=%.18g\n", log_nom);
1846  display_vector(m, m_len, "marginals");
1847  display_vector(x, 2*3*counter, "x");
1848 #endif // DEBUG_FISHER_TABLE
1849 
1850  floatmax_t* log_denom_vec=SG_MALLOC(floatmax_t, counter);
1851  SGVector<floatmax_t>::fill_vector(log_denom_vec, counter, (floatmax_t)0.0);
1852 
1853  for (int32_t k=0; k<counter; k++)
1854  {
1855  for (int32_t j=0; j<3; j++)
1856  {
1857  for (int32_t i=0; i<2; i++)
1858  log_denom_vec[k]+=lgammal(x[i+j*2+k*2*3]+1.0);
1859  }
1860  }
1861 
1862  for (int32_t i=0; i<counter; i++)
1863  log_denom_vec[i]=log_nom-log_denom_vec[i];
1864 
1865 #ifdef DEBUG_FISHER_TABLE
1866  display_vector(log_denom_vec, counter, "log_denom_vec");
1867 #endif // DEBUG_FISHER_TABLE
1868 
1869  float64_t nonrand_p=-CMath::INFTY;
1870  for (int32_t i=0; i<counter; i++)
1871  {
1872  if (log_denom_vec[i]<=prob_table_log)
1873  nonrand_p=CMath::logarithmic_sum(nonrand_p, log_denom_vec[i]);
1874  }
1875 
1876 #ifdef DEBUG_FISHER_TABLE
1877  SG_SPRINT("nonrand_p=%.18g\n", nonrand_p);
1878  SG_SPRINT("exp_nonrand_p=%.18g\n", CMath::exp(nonrand_p));
1879 #endif // DEBUG_FISHER_TABLE
1880  nonrand_p=CMath::exp(nonrand_p);
1881 
1882  SG_FREE(log_denom_vec);
1883  SG_FREE(x);
1884  SG_FREE(m);
1885 
1886  return nonrand_p;
1887 }
1888 
1889 float64_t CStatistics::mutual_info(float64_t* p1, float64_t* p2, int32_t len)
1890 {
1891  double e=0;
1892 
1893  for (int32_t i=0; i<len; i++)
1894  for (int32_t j=0; j<len; j++)
1895  e+=exp(p2[j*len+i])*(p2[j*len+i]-p1[i]-p1[j]);
1896 
1897  return (float64_t)e;
1898 }
1899 
1900 float64_t CStatistics::relative_entropy(float64_t* p, float64_t* q, int32_t len)
1901 {
1902  double e=0;
1903 
1904  for (int32_t i=0; i<len; i++)
1905  e+=exp(p[i])*(p[i]-q[i]);
1906 
1907  return (float64_t)e;
1908 }
1909 
1910 float64_t CStatistics::entropy(float64_t* p, int32_t len)
1911 {
1912  double e=0;
1913 
1914  for (int32_t i=0; i<len; i++)
1915  e-=exp(p[i])*p[i];
1916 
1917  return (float64_t)e;
1918 }
1919 
1920 SGVector<int32_t> CStatistics::sample_indices(int32_t sample_size, int32_t N)
1921 {
1922  REQUIRE(sample_size<N,
1923  "sample size should be less than number of indices\n");
1924  int32_t* idxs=SG_MALLOC(int32_t,N);
1925  int32_t i, rnd;
1926  int32_t* permuted_idxs=SG_MALLOC(int32_t,sample_size);
1927 
1928  // reservoir sampling
1929  for (i=0; i<N; i++)
1930  idxs[i]=i;
1931  for (i=0; i<sample_size; i++)
1932  permuted_idxs[i]=idxs[i];
1933  for (i=sample_size; i<N; i++)
1934  {
1935  rnd=CMath::random(1, i);
1936  if (rnd<sample_size)
1937  permuted_idxs[rnd]=idxs[i];
1938  }
1939  SG_FREE(idxs);
1940 
1941  CMath::qsort(permuted_idxs, sample_size);
1942  return SGVector<int32_t>(permuted_idxs, sample_size);
1943 }
1944 
1945 float64_t CStatistics::dlgamma(float64_t x)
1946 {
1947  x = x+6.0;
1948  float64_t df = 1./(x*x);
1949  df = (((df/240-0.003968253986254)*df+1/120.0)*df-1/120.0)*df;
1950  df = df+log(x)-0.5/x-1.0/(x-1.0)-1.0/(x-2.0)-1.0/
1951  (x-3.0)-1.0/(x-4.0)-1.0/(x-5.0)-1.0/(x-6.0);
1952  return df;
1953 }
1954 
SHOGUN Machine Learning Toolbox - Documentation