PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 1998 Pascal Vincent 00005 // Copyright (C) 1999-2002 Pascal Vincent, Yoshua Bengio and University of Montreal 00006 // 00007 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions are met: 00010 // 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 00014 // 2. Redistributions in binary form must reproduce the above copyright 00015 // notice, this list of conditions and the following disclaimer in the 00016 // documentation and/or other materials provided with the distribution. 00017 // 00018 // 3. The name of the authors may not be used to endorse or promote 00019 // products derived from this software without specific prior written 00020 // permission. 00021 // 00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00032 // 00033 // This file is part of the PLearn library. For more information on the PLearn 00034 // library, go to the PLearn Web site at www.plearn.org 00035 00036 #include <plearn/base/general.h> 00037 #include <plearn/math/pl_erf.h> 00038 // #include <iostream> 00039 #include <cmath> 00040 00041 namespace PLearn { 00042 using namespace std; 00043 00044 #define ITMAX 150 00045 #define EPS 3.0e-7 00046 #define FPMIN 1.0e-30 00047 #define Pi 3.141592653589793 00048 #define Log2Pi 1.837877066409 00049 #define Sqrt2Pi 2.506628274631 00050 00051 static double pl_gammln_cof[7]={ 1.000000000190015 , 00052 76.18009172947146 , 00053 -86.50532032941677 , 00054 24.01409824083091 , 00055 -1.231739572450155 , 00056 0.1208650973866179e-2, 00057 -0.5395239384953e-5 }; 00058 00059 // function gamma returns log(Gamma(z)), where 00060 // Gamma(z) = \int_0^infty t^{z-1}*e^{-t} dt 00061 real pl_gammln(real z) 00062 { 00063 double gz,tmp; 00064 static double gamma = 5.0; 00065 gz = (z+0.5)*pl_log(z+gamma+0.5); 00066 gz -= z+gamma+0.5; 00067 gz += 0.5*Log2Pi; 00068 tmp = pl_gammln_cof[0]; 00069 for(int i=1;i<7;i++) tmp += pl_gammln_cof[i]/(z+i); 00070 gz += pl_log(tmp/z); 00071 return(gz); 00072 } 00073 00074 // returns d(pl_gammln(z))/dz also known as the digamma function 00075 real pl_dgammlndz(real z) 00076 { 00077 real tmp0= pl_gammln_cof[0], 00078 tmp1= 0.0; 00079 for(int i= 1; i<7; ++i) 00080 { 00081 tmp0+= pl_gammln_cof[i]/(z+i); 00082 tmp1-= pl_gammln_cof[i]/((z+i)*(z+i)); 00083 } 00084 return (0.5+z)/(5.5+z)-1 + z*(-tmp0/(z*z) + tmp1/z)/tmp0 + pl_log(5.5+z); 00085 } 00086 00087 00088 // returns the series value of 00089 // the incomplete gamma function 00090 real pl_gser(real a, real x) { 00091 real EPSILON = 1e-7; 00092 real g = pl_gammln(a); 00093 real sum,term; 00094 if (x<0 || a<0) 00095 PLERROR("Error in function pl_gser. Bad argument."); 00096 else if (fast_exact_is_equal(x, 0)) 00097 return 0; 00098 00099 sum = term = 1/a; 00100 for(int i=1;i<ITMAX;i++) { 00101 term *= x/(a+i); 00102 sum += term; 00103 if (term < sum*EPSILON) break; 00104 } 00105 return exp(-x+a*pl_log(x)-g)*sum; 00106 } 00107 00108 00109 // returns the continued fraction representation of 00110 // the incomplete gamma function 00111 real pl_gcf(real a, real x) 00112 { 00113 PLASSERT( !is_missing(a) && !is_missing(x) ); 00114 00115 int i; 00116 real an,b,c,d,del,h; 00117 00118 real gln=pl_gammln(a); 00119 b=x+1.0-a; 00120 c=1.0/FPMIN; 00121 d=1.0/b; 00122 h=d; 00123 for (i=1;i<=ITMAX;i++) { 00124 an = -i*(i-a); 00125 b += 2.0; 00126 d=an*d+b; 00127 if (fabs(d) < FPMIN) d=FPMIN; 00128 c=b+an/c; 00129 if (fabs(c) < FPMIN) c=FPMIN; 00130 d=1.0/d; 00131 del=d*c; 00132 h *= del; 00133 if (fabs(del-1.0) < EPS) break; 00134 } 00135 if (i > ITMAX) { 00136 PLWARNING("\"a\" is too large, ITMAX too small in " 00137 "calling pl_gcf(%f,%f)", a,x); 00138 } 00139 return exp(-x+a*pl_log(x)-(gln))*h; 00140 } 00141 00142 00143 // returns the incomplete gamma function Q(a,x) = 1 - P(a,x) 00144 // it either uses the series or the continued fraction formula 00145 real pl_gammq(real a, real x) { 00146 if (x<0 || a<0) 00147 PLERROR("Error in function gammax. Bad arguments."); 00148 if (x<a+1) 00149 return 1-pl_gser(a,x); 00150 return pl_gcf(a,x); 00151 } 00152 00153 00154 // returns the error function "erf" 00155 real pl_erf(real x) { 00156 #ifdef __GNUC__ 00157 //8.5 time faster in my test then plearn version. 00158 return erf(x); 00159 #else 00160 //it is pl_gcf that take too much time...optimise? 00161 return (x<0?-1:1)*(1-pl_gammq(0.5,x*x)); 00162 #endif 00163 } 00164 00165 00166 //returns the gaussian cumulative function 00167 // For X ~ Normal(0,1), cumulative probability function P(X<x) 00168 real gauss_01_cum(real x) { 00169 return 0.5*(1+pl_erf(x*0.707106781187)); 00170 } 00171 00172 // For X ~ Normal(0,1), inverse of cumulative probability function P(X<x) 00173 // i.e. approximately gauss_01_quantile(gauss_01_cum(x)) ~=~ x 00174 // (the inverse is computed with a binary search, the bisection method) 00175 real gauss_01_quantile(real q) { 00176 #ifdef BOUNDCHECK 00177 if(q<0||q>1) 00178 PLERROR("gauss_01_quantile(q=%f) - q is less then 0 or more then 1",q); 00179 PLASSERT(!is_missing(q)); 00180 #endif 00181 00182 // Handle special cases that can lead to infinite loops below. 00183 if (fast_exact_is_equal(q, real(0))) 00184 return -INFINITY; 00185 else if (fast_exact_is_equal(q, real(1))) 00186 return INFINITY; 00187 00188 // first find a reasonable interval (a,b) s.t. cum(a)<q<cum(b) 00189 real a=-2; 00190 real b=2; 00191 real cum_a=gauss_01_cum(a); 00192 real cum_b=gauss_01_cum(b); 00193 while (cum_a>q) { a*=1.5; cum_a=gauss_01_cum(a); } 00194 while (cum_b<q) { b*=1.5; cum_b=gauss_01_cum(b); } 00195 // then start the bisection loop itself 00196 for (;;) { 00197 real c=0.5*(a+b); 00198 real precision = fabs(b-a); 00199 // PRECISION HERE: 00200 if (precision < 1e-6) 00201 return c; 00202 real cum_c = gauss_01_cum(c); 00203 if (cum_c < q) 00204 a=c; 00205 else 00206 b=c; 00207 } 00208 } 00209 00210 00211 // for X ~ Normal(0,1), return density of X at x 00212 real gauss_01_density(real x) 00213 { 00214 return exp(-0.5*x*x) / Sqrt2Pi; 00215 } 00216 00217 real gauss_01_log_density(real x) 00218 { 00219 return -0.5*x*x - 0.5*Log2Pi; 00220 } 00221 00222 real gauss_log_density_var(real x, real mu, real var) 00223 { 00224 real dx=x-mu; 00225 return -0.5*(dx*dx/var + Log2Pi + pl_log(var)); 00226 00227 } 00228 00229 real gauss_density_var(real x, real mu, real var) { 00230 real dx = x - mu; 00231 return exp(-0.5 * dx * dx / var) / Sqrt2Pi; 00232 } 00233 00234 real gauss_log_density_stddev(real x, real mu, real sigma) 00235 { 00236 real dx = (x-mu) / sigma; 00237 return -0.5*(dx*dx + Log2Pi) - pl_log(sigma); 00238 } 00239 00240 real p_value(real mu, real vn) 00241 { 00242 if (is_missing(mu) || is_missing(vn)) 00243 return MISSING_VALUE; 00244 return 1 - gauss_01_cum(fabs(mu/sqrt(vn))); 00245 } 00246 00247 00248 float gaussQuantiletable[GAUSSQUANTILETABLESIZE]; 00249 00250 PLGaussQuantileInitializer::PLGaussQuantileInitializer() 00251 { 00253 real scaling = 1./(GAUSSQUANTILETABLESIZE-1); 00254 for(int i=0; i<GAUSSQUANTILETABLESIZE; i++) { 00255 gaussQuantiletable[i] = (float) gauss_01_quantile(i*scaling); 00256 } 00257 } 00258 00259 PLGaussQuantileInitializer::~PLGaussQuantileInitializer() {} 00260 00261 PLGaussQuantileInitializer pl_gauss_quantile_initializer; 00262 00263 real fast_gauss_01_quantile(real q) 00264 { 00265 #ifdef BOUNDCHECK 00266 if(q<0||q>1) 00267 PLERROR("fast_gauss_01_quantile(q=%f) - " 00268 "q is less then 0 or more then 1",q); 00269 PLASSERT(!is_missing(q)); 00270 #endif 00271 00272 if(q>0.005&&q<0.995) 00273 { 00274 int i; 00275 DOUBLE_TO_INT( double(q*((GAUSSQUANTILETABLESIZE-1))), i); 00276 return real(gaussQuantiletable[i]); 00277 } 00278 else 00279 return gauss_01_quantile(q); 00280 } 00281 } // end of namespace PLearn 00282 00283 00284 /* 00285 Local Variables: 00286 mode:c++ 00287 c-basic-offset:4 00288 c-file-style:"stroustrup" 00289 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00290 indent-tabs-mode:nil 00291 fill-column:79 00292 End: 00293 */ 00294 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :