35 using namespace Eigen;
38 using namespace shogun;
49 return sum/values.
vlen;
80 result=values[median];
86 if (values[low]>values[high])
87 CMath::CMath::swap(values[low], values[high]);
88 result=values[median];
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]);
100 CMath::swap(values[middle], values[low+1]);
108 while (values[low]>values[l]);
111 while (values[h]>values[low]);
114 CMath::swap(values[l], values[h]);
117 CMath::swap(values[low], values[h]);
145 for (i=1; i<values.
vlen; i++)
160 for (i=0; i<values.
vlen; i++)
165 if (values[i]>maxltguess)
166 maxltguess=values[i];
168 else if (values[i]>guess)
171 if (values[i]<mingtguess)
172 mingtguess=values[i];
177 if (less<=(values.
vlen+1)/2&&greater<=(values.
vlen+1)/2)
179 else if (less>greater)
185 if (less>=(values.
vlen+1)/2)
187 else if (less+equal>=(values.
vlen+1)/2)
197 result=median(copy,
true);
205 bool modify,
bool in_place)
212 return median(as_vector, modify, in_place);
221 float64_t mean=CStatistics::mean(values);
225 sum_squared_diff+=CMath::pow(values.
vector[i]-mean, 2);
227 return sum_squared_diff/(values.
vlen-1);
246 result[j]+=values(i,j);
258 result[j]+=values(j,i);
286 result[j]+=CMath::pow(values(i,j)-mean[j], 2);
298 result[j]+=CMath::pow(values(j,i)-mean[j], 2);
309 return CMath::sqrt(variance(values));
317 var[i]=CMath::sqrt(var[i]);
341 centered,
true,
false, 1.0/(observations.
num_rows-1));
357 int32_t deg=values.
vlen-1;
360 float64_t t=CMath::abs(inverse_student_t(deg, alpha));
363 float64_t std_dev=CStatistics::std_deviation(values);
364 float64_t mean=CStatistics::mean(values);
368 conf_int_low=mean-interval;
369 conf_int_up=mean+interval;
382 if (!(k>0 && greater(p, 0)) && less(p, 1))
384 SG_SERROR(
"CStatistics::inverse_student_t_distribution(): "
388 if (greater(p, 0.25) && less(p, 0.75))
396 z=inverse_incomplete_beta(0.5, 0.5*rk, CMath::abs(z));
397 t=CMath::sqrt(rk*z/(1.0-z));
406 if (greater_equal(p, 0.5))
411 z=inverse_incomplete_beta(0.5*rk, 0.5, 2.0*p);
412 if (less(CMath::MAX_REAL_NUMBER*z, rk))
414 result=rflg*CMath::MAX_REAL_NUMBER;
417 t=CMath::sqrt(rk/z-rk);
449 int32_t breaknewtcycle;
450 int32_t breakihalvecycle;
454 if (!(greater_equal(y, 0) && less_equal(y, 1)))
456 SG_SERROR(
"CStatistics::inverse_incomplete_beta(): "
516 if (less_equal(a, 1.0) || less_equal(b, 1.0))
524 yyy=incomplete_beta(aaa, bbb, x);
532 yp=-inverse_normal_cdf(y);
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));
554 if (less(d, CMath::log(CMath::MIN_REAL_NUMBER)))
559 x=aaa/(aaa+bbb*CMath::exp(d));
560 yyy=incomplete_beta(aaa, bbb, x);
562 if (less(CMath::abs(yp), 0.2))
574 if (mainlooppos==ihalve)
579 mainlooppos=ihalvecycle;
586 if (mainlooppos==ihalvecycle)
595 x=1.0-CMath::MACHINE_EPSILON;
606 yyy=incomplete_beta(aaa, bbb, x);
608 if (less(CMath::abs(yp), dithresh))
614 if (less(CMath::abs(yp), dithresh))
633 di=1.0-(1.0-di)*(1.0-di);
648 if (greater(x0, 0.75))
665 yyy=incomplete_beta(aaa, bbb, x);
677 if (rflg==1 && less(x1, CMath::MACHINE_EPSILON))
709 mainlooppos=ihalvecycle;
714 mainlooppos=breakihalvecycle;
722 if (mainlooppos==breakihalvecycle)
724 if (greater_equal(x0, 1.0))
726 x=1.0-CMath::MACHINE_EPSILON;
729 if (less_equal(x, 0.0))
741 if (mainlooppos==newt)
748 lgm=lgamma(aaa+bbb)-lgamma(aaa)-lgamma(bbb);
750 mainlooppos=newtcycle;
757 if (mainlooppos==newtcycle)
763 yyy=incomplete_beta(aaa, bbb, x);
772 if (greater(yyy, yh))
791 if (equal(x, 1.0) || equal(x, 0.0))
793 mainlooppos=breaknewtcycle;
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)))
801 if (greater(d, CMath::log(CMath::MAX_REAL_NUMBER)))
803 mainlooppos=breaknewtcycle;
809 if (less_equal(xt, x0))
812 xt=x0+0.5*yyy*(x-x0);
813 if (less_equal(xt, 0.0))
815 mainlooppos=breaknewtcycle;
819 if (greater_equal(xt, x1))
822 xt=x1-0.5*yyy*(x1-x);
823 if (greater_equal(xt, 1.0))
825 mainlooppos=breaknewtcycle;
830 if (less(CMath::abs(d/x), 128.0*CMath::MACHINE_EPSILON))
835 mainlooppos=newtcycle;
840 mainlooppos=breaknewtcycle;
848 if (mainlooppos==breaknewtcycle)
850 dithresh=256.0*CMath::MACHINE_EPSILON;
861 if (less_equal(x, CMath::MACHINE_EPSILON))
863 x=1.0-CMath::MACHINE_EPSILON;
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);
894 if (!(greater(a, 0) && greater(b, 0)))
896 SG_SERROR(
"CStatistics::incomplete_beta(): "
900 if (!(greater_equal(x, 0) && less_equal(x, 1)))
902 SG_SERROR(
"CStatistics::incomplete_beta(): "
917 if (less_equal(b*x, 1.0) && less_equal(x, 0.95))
919 result=ibetaf_incompletebetaps(a, b, x, maxgam);
923 if (greater(x, a/(a+b)))
936 if ((flag==1 && less_equal(b*x, 1.0)) && less_equal(x, 0.95))
938 t=ibetaf_incompletebetaps(a, b, x, maxgam);
939 if (less_equal(t, CMath::MACHINE_EPSILON))
941 result=1.0-CMath::MACHINE_EPSILON;
949 y=x*(a+b-2.0)-(a-1.0);
952 w=ibetaf_incompletebetafe(a, b, x, big, biginv);
956 w=ibetaf_incompletebetafe2(a, b, x, big, biginv)/xc;
960 if ((less(a+b, maxgam) && less(CMath::abs(y), maxlog))
961 && less(CMath::abs(t), maxlog))
964 t=t*CMath::pow(x, a);
967 t=t*(tgamma(a+b)/(tgamma(a)*tgamma(b)));
970 if (less_equal(t, CMath::MACHINE_EPSILON))
972 result=1.0-CMath::MACHINE_EPSILON;
985 y=y+t+lgamma(a+b)-lgamma(a)-lgamma(b);
997 if (less_equal(t, CMath::MACHINE_EPSILON))
999 t=1.0-CMath::MACHINE_EPSILON;
1013 return inverse_normal_cdf(y0)*std_dev+mean;
1035 expm2=0.13533528323661269189;
1036 s2pi=2.50662827463100050242;
1037 if (less_equal(y0, 0))
1039 result=-CMath::MAX_REAL_NUMBER;
1042 if (greater_equal(y0, 1))
1044 result=CMath::MAX_REAL_NUMBER;
1049 if (greater(y, 1.0-expm2))
1054 if (greater(y, expm2))
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;
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;
1077 x=CMath::sqrt(-2.0*CMath::log(y));
1078 x0=x-CMath::log(x)/x;
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;
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;
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;
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;
1153 z=CMath::MACHINE_EPSILON*ai;
1154 while (greater(CMath::abs(v), z))
1165 if (less(a+b, maxgam)
1166 && less(CMath::abs(u), CMath::log(CMath::MAX_REAL_NUMBER)))
1168 t=tgamma(a+b)/(tgamma(a)*tgamma(b));
1169 s=s*t*CMath::pow(x, a);
1173 t=lgamma(a+b)-lgamma(a)-lgamma(b)+u+CMath::log(s);
1174 if (less(t, CMath::log(CMath::MIN_REAL_NUMBER)))
1227 thresh=3.0*CMath::MACHINE_EPSILON;
1230 xk=-x*k1*k2/(k3*k4);
1244 if (not_equal(qk, 0))
1248 if (not_equal(r, 0))
1250 t=CMath::abs((ans-r)/r);
1257 if (less(t, thresh))
1269 if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1276 if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1332 thresh=3.0*CMath::MACHINE_EPSILON;
1335 xk=-z*k1*k2/(k3*k4);
1349 if (not_equal(qk, 0))
1353 if (not_equal(r, 0))
1355 t=CMath::abs((ans-r)/r);
1362 if (less(t, thresh))
1374 if (greater(CMath::abs(qk)+CMath::abs(pk), big))
1381 if (less(CMath::abs(qk), biginv) || less(CMath::abs(pk), biginv))
1404 igammaepsilon=0.000000000000001;
1405 if (less_equal(x, 0) || less_equal(a, 0))
1410 if (greater(x, 1) && greater(x, a))
1412 result=1-incomplete_gamma_completed(a, x);
1415 ax=a*CMath::log(x)-x-lgamma(a);
1416 if (less(ax, -709.78271289338399))
1431 while (greater(c/ans, igammaepsilon));
1457 igammaepsilon=0.000000000000001;
1458 igammabignumber=4503599627370496.0;
1459 igammabignumberinv=2.22044604925031308085*0.0000000000000001;
1460 if (less_equal(x, 0) || less_equal(a, 0))
1465 if (less(x, 1) || less(x, a))
1467 result=1-incomplete_gamma(a, x);
1470 ax=a*CMath::log(x)-x-lgamma(a);
1471 if (less(ax, -709.78271289338399))
1493 if (not_equal(qk, 0))
1496 t=CMath::abs((ans-r)/r);
1507 if (greater(CMath::abs(pk), igammabignumber))
1509 pkm2=pkm2*igammabignumberinv;
1510 pkm1=pkm1*igammabignumberinv;
1511 qkm2=qkm2*igammabignumberinv;
1512 qkm1=qkm1*igammabignumberinv;
1515 while (greater(t, igammaepsilon));
1523 return incomplete_gamma(a, x/b);
1531 return inverse_incomplete_gamma_completed(a, 1-p)*b;
1552 igammaepsilon=0.000000000000001;
1553 iinvgammabignumber=4503599627370496.0;
1554 x0=iinvgammabignumber;
1558 dithresh=5*igammaepsilon;
1560 y=1-d-inverse_normal_cdf(y0)*CMath::sqrt(d);
1566 if (greater(x, x0) || less(x, x1))
1571 y=incomplete_gamma_completed(a, x);
1572 if (less(y, yl) || greater(y, yh))
1587 d=(a-1)*CMath::log(x)-x-lgm;
1588 if (less(d, -709.78271289338399))
1595 if (less(CMath::abs(d/x), igammaepsilon))
1603 if (equal(x0, iinvgammabignumber))
1605 if (less_equal(x, 0))
1609 while (equal(x0, iinvgammabignumber))
1612 y=incomplete_gamma_completed(a, x);
1628 y=incomplete_gamma_completed(a, x);
1629 lgm=(x0-x1)/(x1+x0);
1630 if (less(CMath::abs(lgm), dithresh))
1635 if (less(CMath::abs(lgm), dithresh))
1639 if (less_equal(x, 0.0))
1643 if (greater_equal(y, y0))
1695 return 0.5*(error_function_complement(-x/std_dev/1.41421356237309514547));
1699 const float64_t CStatistics::ERFC_CASE1=0.0492;
1701 const float64_t CStatistics::ERFC_CASE2=-11.3137;
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;
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,
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;
1743 else if (x<ERFC_CASE2)
1756 return CMath::log(erfc8_weighted_sum(x))-log_of_2-x*x*0.5;
1760 return CMath::log(normal_cdf(x));
1766 return incomplete_gamma(k/2.0,x/2.0);
1772 return incomplete_beta(d1/2.0,d2/2.0,d1*x/(d1*x+d2));
1779 const float64_t sqrt_of_2=1.41421356237309514547;
1783 0.5641895835477550741253201704,
1784 1.275366644729965952479585264,
1785 5.019049726784267463450058,
1786 6.1602098531096305440906,
1787 7.409740605964741794425,
1788 2.97886562639399288862
1794 2.260528520767326969591866945,
1795 9.396034016235054150430579648,
1796 12.0489519278551290360340491,
1797 17.08144074746600431571095,
1798 9.608965327192787870698,
1799 3.3690752069827527677
1807 num=-x*num/sqrt_of_2+P[i];
1813 den=-x*den/sqrt_of_2+Q[i];
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;
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;
1851 if (greater_equal(x, 10))
1856 result=s*(1-error_function_complement(x));
1868 result=2-error_function_complement(-x);
1873 result=1.0-error_function(x);
1876 if (greater_equal(x, 10))
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;
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;
1910 for (int32_t i=0; i<len; i++)
1913 v.
vector[i]=fishers_exact_test_for_2x3_table(table);
1938 for (int32_t i=0; i<3+2; i++)
1939 log_nom+=lgamma(m[i]+1);
1940 log_nom-=lgamma(n+1.0);
1945 for (int32_t i=0; i<3*2; i++)
1953 int32_t dim1=CMath::min(m[0], m[2]);
1957 for (int32_t k=0; k<=dim1; k++)
1959 for (int32_t l=CMath::max(0.0, m[0]-m[4]-k);
1960 l<=CMath::min(m[0]-k, m[3]); l++)
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];
1974 #ifdef DEBUG_FISHER_TABLE
1977 SG_SPRINT(
"l=%g...%g\n", CMath::max(0.0,m[0]-m[4]-0), CMath::min(m[0]-0, m[3]))
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)
1983 display_vector(m, m_len,
"marginals");
1984 display_vector(x, 2*3*counter,
"x");
1985 #endif // DEBUG_FISHER_TABLE
1990 for (int32_t k=0; k<counter; k++)
1992 for (int32_t j=0; j<3; j++)
1994 for (int32_t i=0; i<2; i++)
1995 log_denom_vec[k]+=lgammal(x[i+j*2+k*2*3]+1.0);
1999 for (int32_t i=0; i<counter; i++)
2000 log_denom_vec[i]=log_nom-log_denom_vec[i];
2002 #ifdef DEBUG_FISHER_TABLE
2003 display_vector(log_denom_vec, counter,
"log_denom_vec");
2004 #endif // DEBUG_FISHER_TABLE
2007 for (int32_t i=0; i<counter; i++)
2009 if (log_denom_vec[i]<=prob_table_log)
2010 nonrand_p=CMath::logarithmic_sum(nonrand_p, log_denom_vec[i]);
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);
2019 SG_FREE(log_denom_vec);
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]);
2041 for (int32_t i=0; i<len; i++)
2042 e+=exp(p[i])*(p[i]-q[i]);
2051 for (int32_t i=0; i<len; i++)
2060 "sample size should be less than number of indices\n");
2061 int32_t* idxs=SG_MALLOC(int32_t,N);
2063 int32_t* permuted_idxs=SG_MALLOC(int32_t,sample_size);
2068 for (i=0; i<sample_size; i++)
2069 permuted_idxs[i]=idxs[i];
2070 for (i=sample_size; i<N; i++)
2072 rnd=CMath::random(1, i);
2073 if (rnd<sample_size)
2074 permuted_idxs[rnd]=idxs[i];
2091 result=CMath::PI/CMath::tan(CMath::PI*x);
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};
2125 result+=coeff[i]*power;
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());
2139 PartialPivLU<MatrixXd> lu(eigen_A);
2140 VectorXd tmp(eigen_A.rows());
2142 for (
index_t idx=0; idx<tmp.rows(); idx++)
2145 VectorXd p=lu.permutationP()*tmp;
2148 for (
index_t idx=0; idx<p.rows(); idx++)
2164 VectorXd u=lu.matrixLU().diagonal();
2167 for (
int idx=0; idx<u.rows(); idx++)
2181 result=u.array().abs().log().sum();
2196 MatrixXd l = llt.matrixL();
2199 VectorXd diag = l.diagonal();
2201 for( int32_t i = 0; i < diag.rows(); ++i ) {
2202 retval += log(diag(i));
2211 typedef SparseMatrix<float64_t> MatrixType;
2214 SimplicialLLT<MatrixType> llt;
2218 MatrixType L=llt.matrixL();
2222 for(
index_t i=0; i<M.rows(); ++i )
2223 retval+=log(L.coeff(i,i));
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");
2238 int32_t dim=mean.
vlen;
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();
2249 MatrixXd U=c.llt().matrixU();
2253 if( precision_matrix )
2257 LDLT<MatrixXd> ldlt;
2264 if( !precision_matrix )
2273 for( int32_t i=0; i<N; ++i )
2283 "CStatistics::sample_from_gaussian(): \
2284 Number of covariance rows must be positive!\n");
2286 "CStatistics::sample_from_gaussian(): \
2287 Number of covariance cols must be positive!\n");
2289 "CStatistics::sample_from_gaussian(): \
2290 Covariance is not initialized!\n");
2292 "CStatistics::sample_from_gaussian(): \
2293 Covariance should be square matrix!\n");
2295 "CStatistics::sample_from_gaussian(): \
2296 Mean and covariance dimension mismatch!\n");
2298 int32_t dim=mean.
vlen;
2301 typedef SparseMatrix<float64_t> MatrixType;
2304 SimplicialLLT<MatrixType> llt;
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();
2316 MatrixType LP=llt.matrixL();
2317 MatrixType UP=llt.matrixU();
2321 if( precision_matrix )
2324 SimplicialLLT<MatrixType> lltUP;
2335 s=llt.permutationPinv()*s;
2340 for( int32_t i=0; i<N; ++i )
2346 #endif //HAVE_EIGEN3
2350 SG_SDEBUG(
"entering CStatistics::fit_sigmoid()\n")
2352 REQUIRE(scores.
vector,
"CStatistics::fit_sigmoid() requires "
2353 "scores vector!\n");
2358 SG_SDEBUG(
"counting number of positive and negative labels\n")
2368 SG_SDEBUG(
"%d pos; %d neg\n", prior1, prior0)
2382 float64_t hiTarget=(prior1+1.0)/(prior1+2.0);
2387 for (
index_t i=0; i<length; ++i)
2398 float64_t b=CMath::log((prior0+1.0)/(prior1+1.0));
2401 for (
index_t i=0; i<length; ++i)
2405 fval+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
2407 fval+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
2413 for (it=0; it<maxiter; ++it)
2415 SG_SDEBUG(
"Iteration %d, a=%f, b=%f, fval=%f\n", it, a, b, fval)
2424 for (
index_t i=0; i<length; ++i)
2431 p=CMath::exp(-fApB)/(1.0+CMath::exp(-fApB));
2432 q=1.0/(1.0+CMath::exp(-fApB));
2436 p=1.0/(1.0+CMath::exp(fApB));
2437 q=CMath::exp(fApB)/(1.0+CMath::exp(fApB));
2441 h11+=scores[i]*scores[i]*d2;
2450 if (CMath::abs(g1)<eps && CMath::abs(g2)<eps)
2462 while (stepsize>=minstep)
2469 for (
index_t i=0; i<length; ++i)
2473 newf+=t[i]*fApB+CMath::log(1+CMath::exp(-fApB));
2475 newf+=(t[i]-1)*fApB+CMath::log(1+CMath::exp(fApB));
2479 if (newf<fval+0.0001*stepsize*gd)
2487 stepsize=stepsize/2.0;
2490 if (stepsize<minstep)
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);
2500 SG_SWARNING(
"CStatistics::fit_sigmoid(): reaching maximal iterations,"
2501 " g1=%f, g2=%f\n", g1, g2);
2504 SG_SDEBUG(
"fitted sigmoid: a=%f, b=%f\n", a, b)
2510 SG_SDEBUG(
"leaving CStatistics::fit_sigmoid()\n")
template class SGSparseMatrix
SGSparseVector< T > * sparse_matrix
array of sparse vectors of size num_vectors
index_t num_vectors
total number of vectors
index_t num_features
total number of features
void remove_column_mean()
This class contains some utilities for Eigen3 Sparse Matrix integration with shogun. Currently it provides a method for converting SGSparseMatrix to Eigen3 SparseMatrix.