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

SHOGUN Machine Learning Toolbox - Documentation