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