SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Statistics.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2011 Heiko Strathmann
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  *
10  * ALGLIB Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier under GPL2+
11  * http://www.alglib.net/
12  * See header file for which functions are taken from ALGLIB (with adjustments
13  * for shogun)
14  */
15 
18 
19 using namespace shogun;
20 
22 {
23  ASSERT(values.vlen>0);
24  ASSERT(values.vector);
25 
26  float64_t sum=0;
27  for (index_t i=0; i<values.vlen; ++i)
28  sum+=values.vector[i];
29 
30  return sum/values.vlen;
31 }
32 
34 {
35  ASSERT(values.vlen>1);
36  ASSERT(values.vector);
37 
39 
40  float64_t sum_squared_diff=0;
41  for (index_t i=0; i<values.vlen; ++i)
42  sum_squared_diff+=CMath::pow(values.vector[i]-mean, 2);
43 
44  return sum_squared_diff/(values.vlen-1);
45 }
46 
48 {
49  return CMath::sqrt(variance(values));
50 }
51 
53  float64_t alpha, float64_t& conf_int_low, float64_t& conf_int_up)
54 {
55  ASSERT(values.vlen>1);
56  ASSERT(values.vector);
57 
58  /* using one sided student t distribution evaluation */
59  alpha=alpha/2;
60 
61  /* degrees of freedom */
62  int32_t deg=values.vlen-1;
63 
64  /* compute absolute value of t-value */
66 
67  /* values for calculating confidence interval */
68  float64_t std_dev=std_deviation(values);
70 
71  /* compute confidence interval */
72  float64_t interval=t*std_dev/CMath::sqrt((float64_t)values.vlen);
73  conf_int_low=mean-interval;
74  conf_int_up=mean+interval;
75 
76  return mean;
77 }
78 
80 {
81  float64_t x;
82  float64_t rk;
83  float64_t z;
84  float64_t f;
85  float64_t tz;
86  float64_t p;
87  float64_t xsqk;
88  int32_t j;
89  float64_t result;
90 
91  ASSERT(k>0);
92  if (t==0)
93  {
94  result=0.5;
95  return result;
96  }
97  if (t<-2.0)
98  {
99  rk=k;
100  z=rk/(rk+t*t);
101  result=0.5*incomplete_beta(0.5*rk, 0.5, z);
102  return result;
103  }
104  if (t<0)
105  {
106  x=-t;
107  }
108  else
109  {
110  x=t;
111  }
112  rk=k;
113  z=1.0+x*x/rk;
114  if (k%2!=0)
115  {
116  xsqk=x/CMath::sqrt(rk);
117  p=CMath::atan(xsqk);
118  if (k>1)
119  {
120  f=1.0;
121  tz=1.0;
122  j=3;
123  while (j<=k-2&&tz/f>CMath::MACHINE_EPSILON)
124  {
125  tz=tz*((j-1)/(z*j));
126  f=f+tz;
127  j=j+2;
128  }
129  p=p+f*xsqk/z;
130  }
131  p=p*2.0/CMath::PI;
132  }
133  else
134  {
135  f=1.0;
136  tz=1.0;
137  j=2;
138  while (j<=k-2&&tz/f>CMath::MACHINE_EPSILON)
139  {
140  tz=tz*((j-1)/(z*j));
141  f=f+tz;
142  j=j+2;
143  }
144  p=f*x/CMath::sqrt(z*rk);
145  }
146  if (t<0)
147  {
148  p=-p;
149  }
150  result=0.5+0.5*p;
151  return result;
152 }
153 
155 {
156  float64_t t;
157  float64_t xc;
158  float64_t w;
159  float64_t y;
160  int32_t flag;
161  float64_t big;
162  float64_t biginv;
163  float64_t maxgam;
164  float64_t minlog;
165  float64_t maxlog;
166  float64_t result;
167 
168  big=4.503599627370496e15;
169  biginv=2.22044604925031308085e-16;
170  maxgam=171.624376956302725;
173  ASSERT(a>0&&b>0);
174  ASSERT(x>=0&&x<=1);
175  if (x==0)
176  {
177  result=0;
178  return result;
179  }
180  if (x==1)
181  {
182  result=1;
183  return result;
184  }
185  flag=0;
186  if (b*x<=1.0&&x<=0.95)
187  {
188  result=ibetaf_incomplete_beta_ps(a, b, x, maxgam);
189  return result;
190  }
191  w=1.0-x;
192  if (x>a/(a+b))
193  {
194  flag=1;
195  t=a;
196  a=b;
197  b=t;
198  xc=x;
199  x=w;
200  }
201  else
202  {
203  xc=w;
204  }
205  if ((flag==1&&b*x<=1.0)&&x<=0.95)
206  {
207  t=ibetaf_incomplete_beta_ps(a, b, x, maxgam);
208  if (t<=CMath::MACHINE_EPSILON)
209  {
210  result=1.0-CMath::MACHINE_EPSILON;
211  }
212  else
213  {
214  result=1.0-t;
215  }
216  return result;
217  }
218  y=x*(a+b-2.0)-(a-1.0);
219  if (y<0.0)
220  {
221  w=ibetaf_incomplete_beta_fe(a, b, x, big, biginv);
222  }
223  else
224  {
225  w=ibetaf_incomplete_beta_fe2(a, b, x, big, biginv)/xc;
226  }
227  y=a*CMath::log(x);
228  t=b*CMath::log(xc);
229  if ((a+b<maxgam&&CMath::abs(y)<maxlog)&&CMath::abs(t)<maxlog)
230  {
231  t=CMath::pow(xc, b);
232  t=t*CMath::pow(x, a);
233  t=t/a;
234  t=t*w;
235  t=t*(CMath::tgamma(a+b)/(CMath::tgamma(a)*CMath::tgamma(b)));
236  if (flag==1)
237  {
238  if (t<=CMath::MACHINE_EPSILON)
239  {
240  result=1.0-CMath::MACHINE_EPSILON;
241  }
242  else
243  {
244  result=1.0-t;
245  }
246  }
247  else
248  {
249  result=t;
250  }
251  return result;
252  }
253  y=y+t+CMath::lgamma(a+b)-CMath::lgamma(a)-CMath::lgamma(b);
254  y=y+CMath::log(w/a);
255  if (y<minlog)
256  {
257  t=0.0;
258  }
259  else
260  {
261  t=CMath::exp(y);
262  }
263  if (flag==1)
264  {
265  if (t<=CMath::MACHINE_EPSILON)
266  {
268  }
269  else
270  {
271  t=1.0-t;
272  }
273  }
274  result=t;
275  return result;
276 }
277 
279  float64_t x, float64_t maxgam)
280 {
281  float64_t s;
282  float64_t t;
283  float64_t u;
284  float64_t v;
285  float64_t n;
286  float64_t t1;
287  float64_t z;
288  float64_t ai;
289  float64_t result;
290 
291  ai=1.0/a;
292  u=(1.0-b)*x;
293  v=u/(a+1.0);
294  t1=v;
295  t=u;
296  n=2.0;
297  s=0.0;
299  while (CMath::abs(v)>z)
300  {
301  u=(n-b)*x/n;
302  t=t*u;
303  v=t/(a+n);
304  s=s+v;
305  n=n+1.0;
306  }
307  s=s+t1;
308  s=s+ai;
309  u=a*CMath::log(x);
310  if (a+b<maxgam&&CMath::abs(u)<CMath::log(CMath::MAX_REAL_NUMBER))
311  {
313  s=s*t*CMath::pow(x, a);
314  }
315  else
316  {
317  t=CMath::lgamma(a+b)-CMath::lgamma(a)-CMath::lgamma(b)+u+CMath::log(s);
318  if (t<CMath::log(CMath::MIN_REAL_NUMBER))
319  {
320  s=0.0;
321  }
322  else
323  {
324  s=CMath::exp(t);
325  }
326  }
327  result=s;
328  return result;
329 }
330 
332  float64_t x, float64_t big, float64_t biginv)
333 {
334  float64_t xk;
335  float64_t pk;
336  float64_t pkm1;
337  float64_t pkm2;
338  float64_t qk;
339  float64_t qkm1;
340  float64_t qkm2;
341  float64_t k1;
342  float64_t k2;
343  float64_t k3;
344  float64_t k4;
345  float64_t k5;
346  float64_t k6;
347  float64_t k7;
348  float64_t k8;
349  float64_t r;
350  float64_t t;
351  float64_t ans;
352  float64_t z;
353  float64_t thresh;
354  int32_t n;
355  float64_t result;
356 
357  k1=a;
358  k2=b-1.0;
359  k3=a;
360  k4=a+1.0;
361  k5=1.0;
362  k6=a+b;
363  k7=a+1.0;
364  k8=a+2.0;
365  pkm2=0.0;
366  qkm2=1.0;
367  pkm1=1.0;
368  qkm1=1.0;
369  z=x/(1.0-x);
370  ans=1.0;
371  r=1.0;
372  n=0;
373  thresh=3.0*CMath::MACHINE_EPSILON;
374  do
375  {
376  xk=-z*k1*k2/(k3*k4);
377  pk=pkm1+pkm2*xk;
378  qk=qkm1+qkm2*xk;
379  pkm2=pkm1;
380  pkm1=pk;
381  qkm2=qkm1;
382  qkm1=qk;
383  xk=z*k5*k6/(k7*k8);
384  pk=pkm1+pkm2*xk;
385  qk=qkm1+qkm2*xk;
386  pkm2=pkm1;
387  pkm1=pk;
388  qkm2=qkm1;
389  qkm1=qk;
390  if (qk!=0)
391  {
392  r=pk/qk;
393  }
394  if (r!=0)
395  {
396  t=CMath::abs((ans-r)/r);
397  ans=r;
398  }
399  else
400  {
401  t=1.0;
402  }
403  if (t<thresh)
404  {
405  break;
406  }
407  k1=k1+1.0;
408  k2=k2-1.0;
409  k3=k3+2.0;
410  k4=k4+2.0;
411  k5=k5+1.0;
412  k6=k6+1.0;
413  k7=k7+2.0;
414  k8=k8+2.0;
415  if (CMath::abs(qk)+CMath::abs(pk)>big)
416  {
417  pkm2=pkm2*biginv;
418  pkm1=pkm1*biginv;
419  qkm2=qkm2*biginv;
420  qkm1=qkm1*biginv;
421  }
422  if (CMath::abs(qk)<biginv||CMath::abs(pk)<biginv)
423  {
424  pkm2=pkm2*big;
425  pkm1=pkm1*big;
426  qkm2=qkm2*big;
427  qkm1=qkm1*big;
428  }
429  n=n+1;
430  } while (n!=300);
431  result=ans;
432  return result;
433 }
434 
436  float64_t x, float64_t big, float64_t biginv)
437 {
438  float64_t xk;
439  float64_t pk;
440  float64_t pkm1;
441  float64_t pkm2;
442  float64_t qk;
443  float64_t qkm1;
444  float64_t qkm2;
445  float64_t k1;
446  float64_t k2;
447  float64_t k3;
448  float64_t k4;
449  float64_t k5;
450  float64_t k6;
451  float64_t k7;
452  float64_t k8;
453  float64_t r;
454  float64_t t;
455  float64_t ans;
456  float64_t thresh;
457  int32_t n;
458  float64_t result;
459 
460  k1=a;
461  k2=a+b;
462  k3=a;
463  k4=a+1.0;
464  k5=1.0;
465  k6=b-1.0;
466  k7=k4;
467  k8=a+2.0;
468  pkm2=0.0;
469  qkm2=1.0;
470  pkm1=1.0;
471  qkm1=1.0;
472  ans=1.0;
473  r=1.0;
474  n=0;
475  thresh=3.0*CMath::MACHINE_EPSILON;
476  do
477  {
478  xk=-x*k1*k2/(k3*k4);
479  pk=pkm1+pkm2*xk;
480  qk=qkm1+qkm2*xk;
481  pkm2=pkm1;
482  pkm1=pk;
483  qkm2=qkm1;
484  qkm1=qk;
485  xk=x*k5*k6/(k7*k8);
486  pk=pkm1+pkm2*xk;
487  qk=qkm1+qkm2*xk;
488  pkm2=pkm1;
489  pkm1=pk;
490  qkm2=qkm1;
491  qkm1=qk;
492  if (qk!=0)
493  {
494  r=pk/qk;
495  }
496  if (r!=0)
497  {
498  t=CMath::abs((ans-r)/r);
499  ans=r;
500  }
501  else
502  {
503  t=1.0;
504  }
505  if (t<thresh)
506  {
507  break;
508  }
509  k1=k1+1.0;
510  k2=k2+1.0;
511  k3=k3+2.0;
512  k4=k4+2.0;
513  k5=k5+1.0;
514  k6=k6-1.0;
515  k7=k7+2.0;
516  k8=k8+2.0;
517  if (CMath::abs(qk)+CMath::abs(pk)>big)
518  {
519  pkm2=pkm2*biginv;
520  pkm1=pkm1*biginv;
521  qkm2=qkm2*biginv;
522  qkm1=qkm1*biginv;
523  }
524  if (CMath::abs(qk)<biginv||CMath::abs(pk)<biginv)
525  {
526  pkm2=pkm2*big;
527  pkm1=pkm1*big;
528  qkm2=qkm2*big;
529  qkm1=qkm1*big;
530  }
531  n=n+1;
532  } while (n!=300);
533  result=ans;
534  return result;
535 }
536 
538 {
539  float64_t t;
540  float64_t rk;
541  float64_t z;
542  int32_t rflg;
543  float64_t result;
544 
545  ASSERT(k>0&&p>0&&p<1);
546  rk=k;
547  if (p>0.25&&p<0.75)
548  {
549  if (p==0.5)
550  {
551  result=0;
552  return result;
553  }
554  z=1.0-2.0*p;
555  z=inverse_incomplete_beta(0.5, 0.5*rk, CMath::abs(z));
556  t=CMath::sqrt(rk*z/(1.0-z));
557  if (p<0.5)
558  {
559  t=-t;
560  }
561  result=t;
562  return result;
563  }
564  rflg=-1;
565  if (p>=0.5)
566  {
567  p=1.0-p;
568  rflg=1;
569  }
570  z=inverse_incomplete_beta(0.5*rk, 0.5, 2.0*p);
571  if (CMath::MAX_REAL_NUMBER*z<rk)
572  {
573  result=rflg*CMath::MAX_REAL_NUMBER;
574  return result;
575  }
576  t=CMath::sqrt(rk/z-rk);
577  result=rflg*t;
578  return result;
579 }
580 
582  float64_t y)
583 {
584  float64_t aaa;
585  float64_t bbb;
586  float64_t y0;
587  float64_t d;
588  float64_t yyy;
589  float64_t x;
590  float64_t x0;
591  float64_t x1;
592  float64_t lgm;
593  float64_t yp;
594  float64_t di;
595  float64_t dithresh;
596  float64_t yl;
597  float64_t yh;
598  float64_t xt;
599  int32_t i;
600  int32_t rflg;
601  int32_t dir;
602  int32_t nflg;
603  int32_t mainlooppos;
604  int32_t ihalve;
605  int32_t ihalvecycle;
606  int32_t newt;
607  int32_t newtcycle;
608  int32_t breaknewtcycle;
609  int32_t breakihalvecycle;
610  float64_t result;
611 
612  i=0;
613  ASSERT(y>=0&&y<=1);
614 
615  /*
616  * special cases
617  */
618  if (y==0)
619  {
620  result=0;
621  return result;
622  }
623  if (y==1.0)
624  {
625  result=1;
626  return result;
627  }
628 
629  /*
630  * these initializations are not really necessary,
631  * but without them compiler complains about 'possibly uninitialized variables'.
632  */
633  dithresh=0;
634  rflg=0;
635  aaa=0;
636  bbb=0;
637  y0=0;
638  x=0;
639  yyy=0;
640  lgm=0;
641  dir=0;
642  di=0;
643 
644  /*
645  * normal initializations
646  */
647  x0=0.0;
648  yl=0.0;
649  x1=1.0;
650  yh=1.0;
651  nflg=0;
652  mainlooppos=0;
653  ihalve=1;
654  ihalvecycle=2;
655  newt=3;
656  newtcycle=4;
657  breaknewtcycle=5;
658  breakihalvecycle=6;
659 
660  /*
661  * main loop
662  */
663  for (;;)
664  {
665 
666  /*
667  * start
668  */
669  if (mainlooppos==0)
670  {
671  if (a<=1.0||b<=1.0)
672  {
673  dithresh=1.0e-6;
674  rflg=0;
675  aaa=a;
676  bbb=b;
677  y0=y;
678  x=aaa/(aaa+bbb);
679  yyy=incomplete_beta(aaa, bbb, x);
680  mainlooppos=ihalve;
681  continue;
682  }
683  else
684  {
685  dithresh=1.0e-4;
686  }
688  if (y>0.5)
689  {
690  rflg=1;
691  aaa=b;
692  bbb=a;
693  y0=1.0-y;
694  yp=-yp;
695  }
696  else
697  {
698  rflg=0;
699  aaa=a;
700  bbb=b;
701  y0=y;
702  }
703  lgm=(yp*yp-3.0)/6.0;
704  x=2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
705  d=yp*CMath::sqrt(x+lgm)/x
706  -(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))
707  *(lgm+5.0/6.0-2.0/(3.0*x));
708  d=2.0*d;
710  {
711  x=0;
712  break;
713  }
714  x=aaa/(aaa+bbb*CMath::exp(d));
715  yyy=incomplete_beta(aaa, bbb, x);
716  yp=(yyy-y0)/y0;
717  if (CMath::abs(yp)<0.2)
718  {
719  mainlooppos=newt;
720  continue;
721  }
722  mainlooppos=ihalve;
723  continue;
724  }
725 
726  /*
727  * ihalve
728  */
729  if (mainlooppos==ihalve)
730  {
731  dir=0;
732  di=0.5;
733  i=0;
734  mainlooppos=ihalvecycle;
735  continue;
736  }
737 
738  /*
739  * ihalvecycle
740  */
741  if (mainlooppos==ihalvecycle)
742  {
743  if (i<=99)
744  {
745  if (i!=0)
746  {
747  x=x0+di*(x1-x0);
748  if (x==1.0)
749  {
751  }
752  if (x==0.0)
753  {
754  di=0.5;
755  x=x0+di*(x1-x0);
756  if (x==0.0)
757  {
758  break;
759  }
760  }
761  yyy=incomplete_beta(aaa, bbb, x);
762  yp=(x1-x0)/(x1+x0);
763  if (CMath::abs(yp)<dithresh)
764  {
765  mainlooppos=newt;
766  continue;
767  }
768  yp=(yyy-y0)/y0;
769  if (CMath::abs(yp)<dithresh)
770  {
771  mainlooppos=newt;
772  continue;
773  }
774  }
775  if (yyy<y0)
776  {
777  x0=x;
778  yl=yyy;
779  if (dir<0)
780  {
781  dir=0;
782  di=0.5;
783  }
784  else
785  {
786  if (dir>3)
787  {
788  di=1.0-(1.0-di)*(1.0-di);
789  }
790  else
791  {
792  if (dir>1)
793  {
794  di=0.5*di+0.5;
795  }
796  else
797  {
798  di=(y0-yyy)/(yh-yl);
799  }
800  }
801  }
802  dir=dir+1;
803  if (x0>0.75)
804  {
805  if (rflg==1)
806  {
807  rflg=0;
808  aaa=a;
809  bbb=b;
810  y0=y;
811  }
812  else
813  {
814  rflg=1;
815  aaa=b;
816  bbb=a;
817  y0=1.0-y;
818  }
819  x=1.0-x;
820  yyy=incomplete_beta(aaa, bbb, x);
821  x0=0.0;
822  yl=0.0;
823  x1=1.0;
824  yh=1.0;
825  mainlooppos=ihalve;
826  continue;
827  }
828  }
829  else
830  {
831  x1=x;
832  if (rflg==1&&x1<CMath::MACHINE_EPSILON)
833  {
834  x=0.0;
835  break;
836  }
837  yh=yyy;
838  if (dir>0)
839  {
840  dir=0;
841  di=0.5;
842  }
843  else
844  {
845  if (dir<-3)
846  {
847  di=di*di;
848  }
849  else
850  {
851  if (dir<-1)
852  {
853  di=0.5*di;
854  }
855  else
856  {
857  di=(yyy-y0)/(yh-yl);
858  }
859  }
860  }
861  dir=dir-1;
862  }
863  i=i+1;
864  mainlooppos=ihalvecycle;
865  continue;
866  }
867  else
868  {
869  mainlooppos=breakihalvecycle;
870  continue;
871  }
872  }
873 
874  /*
875  * breakihalvecycle
876  */
877  if (mainlooppos==breakihalvecycle)
878  {
879  if (x0>=1.0)
880  {
882  break;
883  }
884  if (x<=0.0)
885  {
886  x=0.0;
887  break;
888  }
889  mainlooppos=newt;
890  continue;
891  }
892 
893  /*
894  * newt
895  */
896  if (mainlooppos==newt)
897  {
898  if (nflg!=0)
899  {
900  break;
901  }
902  nflg=1;
903  lgm=CMath::lgamma(aaa+bbb)-CMath::lgamma(aaa)-CMath::lgamma(bbb);
904  i=0;
905  mainlooppos=newtcycle;
906  continue;
907  }
908 
909  /*
910  * newtcycle
911  */
912  if (mainlooppos==newtcycle)
913  {
914  if (i<=7)
915  {
916  if (i!=0)
917  {
918  yyy=incomplete_beta(aaa, bbb, x);
919  }
920  if (yyy<yl)
921  {
922  x=x0;
923  yyy=yl;
924  }
925  else
926  {
927  if (yyy>yh)
928  {
929  x=x1;
930  yyy=yh;
931  }
932  else
933  {
934  if (yyy<y0)
935  {
936  x0=x;
937  yl=yyy;
938  }
939  else
940  {
941  x1=x;
942  yh=yyy;
943  }
944  }
945  }
946  if (x==1.0||x==0.0)
947  {
948  mainlooppos=breaknewtcycle;
949  continue;
950  }
951  d=(aaa-1.0)*CMath::log(x)+(bbb-1.0)*CMath::log(1.0-x)+lgm;
953  {
954  break;
955  }
957  {
958  mainlooppos=breaknewtcycle;
959  continue;
960  }
961  d=CMath::exp(d);
962  d=(yyy-y0)/d;
963  xt=x-d;
964  if (xt<=x0)
965  {
966  yyy=(x-x0)/(x1-x0);
967  xt=x0+0.5*yyy*(x-x0);
968  if (xt<=0.0)
969  {
970  mainlooppos=breaknewtcycle;
971  continue;
972  }
973  }
974  if (xt>=x1)
975  {
976  yyy=(x1-x)/(x1-x0);
977  xt=x1-0.5*yyy*(x1-x);
978  if (xt>=1.0)
979  {
980  mainlooppos=breaknewtcycle;
981  continue;
982  }
983  }
984  x=xt;
985  if (CMath::abs(d/x)<128.0*CMath::MACHINE_EPSILON)
986  {
987  break;
988  }
989  i=i+1;
990  mainlooppos=newtcycle;
991  continue;
992  }
993  else
994  {
995  mainlooppos=breaknewtcycle;
996  continue;
997  }
998  }
999 
1000  /*
1001  * breaknewtcycle
1002  */
1003  if (mainlooppos==breaknewtcycle)
1004  {
1005  dithresh=256.0*CMath::MACHINE_EPSILON;
1006  mainlooppos=ihalve;
1007  continue;
1008  }
1009  }
1010 
1011  /*
1012  * done
1013  */
1014  if (rflg!=0)
1015  {
1016  if (x<=CMath::MACHINE_EPSILON)
1017  {
1018  x=1.0-CMath::MACHINE_EPSILON;
1019  }
1020  else
1021  {
1022  x=1.0-x;
1023  }
1024  }
1025  result=x;
1026  return result;
1027 }
1028 
1030 {
1031  float64_t expm2;
1032  float64_t s2pi;
1033  float64_t x;
1034  float64_t y;
1035  float64_t z;
1036  float64_t y2;
1037  float64_t x0;
1038  float64_t x1;
1039  int32_t code;
1040  float64_t p0;
1041  float64_t q0;
1042  float64_t p1;
1043  float64_t q1;
1044  float64_t p2;
1045  float64_t q2;
1046  float64_t result;
1047 
1048  expm2=0.13533528323661269189;
1049  s2pi=2.50662827463100050242;
1050  if (y0<=0)
1051  {
1052  result=-CMath::MAX_REAL_NUMBER;
1053  return result;
1054  }
1055  if (y0>=1)
1056  {
1057  result=CMath::MAX_REAL_NUMBER;
1058  return result;
1059  }
1060  code=1;
1061  y=y0;
1062  if (y>1.0-expm2)
1063  {
1064  y=1.0-y;
1065  code=0;
1066  }
1067  if (y>expm2)
1068  {
1069  y=y-0.5;
1070  y2=y*y;
1071  p0=-59.9633501014107895267;
1072  p0=98.0010754185999661536+y2*p0;
1073  p0=-56.6762857469070293439+y2*p0;
1074  p0=13.9312609387279679503+y2*p0;
1075  p0=-1.23916583867381258016+y2*p0;
1076  q0=1;
1077  q0=1.95448858338141759834+y2*q0;
1078  q0=4.67627912898881538453+y2*q0;
1079  q0=86.3602421390890590575+y2*q0;
1080  q0=-225.462687854119370527+y2*q0;
1081  q0=200.260212380060660359+y2*q0;
1082  q0=-82.0372256168333339912+y2*q0;
1083  q0=15.9056225126211695515+y2*q0;
1084  q0=-1.18331621121330003142+y2*q0;
1085  x=y+y*y2*p0/q0;
1086  x=x*s2pi;
1087  result=x;
1088  return result;
1089  }
1090  x=CMath::sqrt(-2.0*CMath::log(y));
1091  x0=x-CMath::log(x)/x;
1092  z=1.0/x;
1093  if (x<8.0)
1094  {
1095  p1=4.05544892305962419923;
1096  p1=31.5251094599893866154+z*p1;
1097  p1=57.1628192246421288162+z*p1;
1098  p1=44.0805073893200834700+z*p1;
1099  p1=14.6849561928858024014+z*p1;
1100  p1=2.18663306850790267539+z*p1;
1101  p1=-1.40256079171354495875*0.1+z*p1;
1102  p1=-3.50424626827848203418*0.01+z*p1;
1103  p1=-8.57456785154685413611*0.0001+z*p1;
1104  q1=1;
1105  q1=15.7799883256466749731+z*q1;
1106  q1=45.3907635128879210584+z*q1;
1107  q1=41.3172038254672030440+z*q1;
1108  q1=15.0425385692907503408+z*q1;
1109  q1=2.50464946208309415979+z*q1;
1110  q1=-1.42182922854787788574*0.1+z*q1;
1111  q1=-3.80806407691578277194*0.01+z*q1;
1112  q1=-9.33259480895457427372*0.0001+z*q1;
1113  x1=z*p1/q1;
1114  }
1115  else
1116  {
1117  p2=3.23774891776946035970;
1118  p2=6.91522889068984211695+z*p2;
1119  p2=3.93881025292474443415+z*p2;
1120  p2=1.33303460815807542389+z*p2;
1121  p2=2.01485389549179081538*0.1+z*p2;
1122  p2=1.23716634817820021358*0.01+z*p2;
1123  p2=3.01581553508235416007*0.0001+z*p2;
1124  p2=2.65806974686737550832*0.000001+z*p2;
1125  p2=6.23974539184983293730*0.000000001+z*p2;
1126  q2=1;
1127  q2=6.02427039364742014255+z*q2;
1128  q2=3.67983563856160859403+z*q2;
1129  q2=1.37702099489081330271+z*q2;
1130  q2=2.16236993594496635890*0.1+z*q2;
1131  q2=1.34204006088543189037*0.01+z*q2;
1132  q2=3.28014464682127739104*0.0001+z*q2;
1133  q2=2.89247864745380683936*0.000001+z*q2;
1134  q2=6.79019408009981274425*0.000000001+z*q2;
1135  x1=z*p2/q2;
1136  }
1137  x=x0-x1;
1138  if (code!=0)
1139  {
1140  x=-x;
1141  }
1142  result=x;
1143  return result;
1144 }

SHOGUN Machine Learning Toolbox - Documentation