PLearn 0.1
pl_erf.cc
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines