Autor Tema: Integración de la Campana de Gauss (Leído 8721 veces)

PicMinor · « **en:** 26 de Mayo de 2014, 06:07:35 »

¡ Saludos al foro !

En estadística se utiliza la integral de la campana de Gauss para calcular probabilidades. En los libros siempre hablan de usar las Tablas pero para hacer los cálculos por ordenador no he localizado ningún algoritmo que genere los valores de las tablas. ¿Alguien sabe de alguno?

¡ Gracias de antemano!

Picuino · « **Respuesta #1 en:** 26 de Mayo de 2014, 07:05:27 »

Creo que lo que buscas es la función error.
http://es.wikipedia.org/wiki/Funci%C3%B3n_error

Lo tengo en algún directorio del disco duro.
Te lo busco y lo posteo.

Saludos.

Picuino · « **Respuesta #2 en:** 26 de Mayo de 2014, 07:29:31 »

Unos enlaces interesantes:

http://www.thefreecountry.com/sourcecode/mathematics.shtml

http://www.jhauser.us/arithmetic/SoftFloat.html

http://www.gnu.org/software/libc/index.html

http://netlib.org/cephes/

Saludos.

Picuino · « **Respuesta #3 en:** 26 de Mayo de 2014, 08:49:50 »

FUNCION ERF.C

http://www.ks.uiuc.edu/Research/namd/doxygen/erf_8C-source.html
http://www.johndcook.com/cpp_erf.html

Código: [Seleccionar]

#include <cmath>

double erf(double x)
{
    // constants
    double a1 =  0.254829592;
    double a2 = -0.284496736;
    double a3 =  1.421413741;
    double a4 = -1.453152027;
    double a5 =  1.061405429;
    double p  =  0.3275911;

    // Save the sign of x
    int sign = 1;
    if (x < 0)
        sign = -1;
    x = fabs(x);

    // A&S formula 7.1.26
    double t = 1.0/(1.0 + p*x);
    double y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return sign*y;
}

La fórmula está tomada de "Abramowitz and Stegun: Handbook of Mathematical Functions" 7.1.26
http://people.math.sfu.ca/~cbm/aands/page_299.htm

Saludos.

PicMinor · « **Respuesta #4 en:** 26 de Mayo de 2014, 11:21:22 »

¡Saludos al Foro!

Para Picuino: El código que me has puesto es de la función error. Esta función toma valores entre -1 y 1 y pasa por 0,0. La integral de la distribución normal va de 0 a 1 y pasa por (0,0.5). El código que has puesto es similar a la aproximación de Hastings que se puede ver en este enlace:

Aproximación de Hastings

Pero verás que las constantes que utiliza son del tipo Single. Estoy buscando algún algoritmo recursivo o polinomial que me permita obtener precisión del orden del Double.

¡Gracias de antemano!

Picuino · « **Respuesta #5 en:** 26 de Mayo de 2014, 16:54:31 »

La Campana de Gauss o función gaussiana es esta:

La fórmula que la genera es esta:

En el caso de que la constante a valga:

a =

la función gaussiana será la función de densidad de la distribución normal.

La Distribución Normal es la integral de la campana de Gauss, en el caso que comenté de que:

a =

La Distribución Normal tiene este aspecto:

La fórmula de la Distribución Normal es:

Ambas fórmulas están definidas para x en el rango de menos infinito hasta más infinito.
En las librerías matemáticas no es usual encontrar la Distribución Normal, en cambio si suele estar presente la función de error.

La función de error tiene este aspecto:

Y su fórmula es:

Como puedes ver, no es exactamente igual a la Distribución Normal. Pero no es difícil relacionarlas.
Como explica en el enlace que te dejé, en el apartado de funciones relacionadas:

http://es.wikipedia.org/wiki/Funci%C3%B3n_error#Funciones_relacionadas
"La función error es esencialmente idéntica a la función distribución de probabilidad normal estándar, designada como Φ, ya que su única diferencia es su escala y una traslación."

La mayoría de los lenguajes serios tendrá una biblioteca en la que aparezca la función de error con muchos decimales de exactitud. Si no tienes la biblioteca, puedes usar la fórmula que comentabas (Aproximación de Hastings) o la que yo te dejé.

¿Para qué quieres más de 6 dígitos de exactitud?

¿Por qué no utilizar la función error que ya viene con el lenguaje?

Saludos.

Picuino · « **Respuesta #6 en:** 26 de Mayo de 2014, 17:07:11 »

Buscando en Google "erf source code double"

Sale esto:

Código: [Seleccionar]

00001 /*
00002  * Copied from OpenBSD project (src/lib/libm/src/s_erf.c)
00003  * Specialized for 32-bit little endian architectures.
00004  */
00005 
00006 /* Real math libraries provide erf(), CUDA also provides an implementation. */
00007 #if defined(WIN32) && !defined(NAMD_CUDA)
00008 
00009 /*
00010  * ====================================================
00011  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00012  *
00013  * Developed at SunPro, a Sun Microsystems, Inc. business.
00014  * Permission to use, copy, modify, and distribute this
00015  * software is freely granted, provided that this notice 
00016  * is preserved.
00017  * ====================================================
00018  */
00019 
00020 /* double erf(double x)
00021  * double erfc(double x)
00022  *                           x
00023  *                    2      |\
00024  *     erf(x)  =  ---------  | exp(-t*t)dt
00025  *                 sqrt(pi) \| 
00026  *                           0
00027  *
00028  *     erfc(x) =  1-erf(x)
00029  *  Note that 
00030  *              erf(-x) = -erf(x)
00031  *              erfc(-x) = 2 - erfc(x)
00032  *
00033  * Method:
00034  *      1. For |x| in [0, 0.84375]
00035  *          erf(x)  = x + x*R(x^2)
00036  *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
00037  *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
00038  *         where R = P/Q where P is an odd poly of degree 8 and
00039  *         Q is an odd poly of degree 10.
00040  *                                               -57.90
00041  *                      | R - (erf(x)-x)/x | <= 2
00042  *      
00043  *
00044  *         Remark. The formula is derived by noting
00045  *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
00046  *         and that
00047  *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
00048  *         is close to one. The interval is chosen because the fix
00049  *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
00050  *         near 0.6174), and by some experiment, 0.84375 is chosen to
00051  *         guarantee the error is less than one ulp for erf.
00052  *
00053  *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
00054  *         c = 0.84506291151 rounded to single (24 bits)
00055  *              erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
00056  *              erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
00057  *                        1+(c+P1(s)/Q1(s))    if x < 0
00058  *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
00059  *         Remark: here we use the taylor series expansion at x=1.
00060  *              erf(1+s) = erf(1) + s*Poly(s)
00061  *                       = 0.845.. + P1(s)/Q1(s)
00062  *         That is, we use rational approximation to approximate
00063  *                      erf(1+s) - (c = (single)0.84506291151)
00064  *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
00065  *         where 
00066  *              P1(s) = degree 6 poly in s
00067  *              Q1(s) = degree 6 poly in s
00068  *
00069  *      3. For x in [1.25,1/0.35(~2.857143)], 
00070  *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
00071  *              erf(x)  = 1 - erfc(x)
00072  *         where 
00073  *              R1(z) = degree 7 poly in z, (z=1/x^2)
00074  *              S1(z) = degree 8 poly in z
00075  *
00076  *      4. For x in [1/0.35,28]
00077  *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
00078  *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
00079  *                      = 2.0 - tiny            (if x <= -6)
00080  *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
00081  *              erf(x)  = sign(x)*(1.0 - tiny)
00082  *         where
00083  *              R2(z) = degree 6 poly in z, (z=1/x^2)
00084  *              S2(z) = degree 7 poly in z
00085  *
00086  *      Note1:
00087  *         To compute exp(-x*x-0.5625+R/S), let s be a single
00088  *         precision number and s := x; then
00089  *              -x*x = -s*s + (s-x)*(s+x)
00090  *              exp(-x*x-0.5626+R/S) = 
00091  *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
00092  *      Note2:
00093  *         Here 4 and 5 make use of the asymptotic series
00094  *                        exp(-x*x)
00095  *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
00096  *                        x*sqrt(pi)
00097  *         We use rational approximation to approximate
00098  *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
00099  *         Here is the error bound for R1/S1 and R2/S2
00100  *              |R1/S1 - f(x)|  < 2**(-62.57)
00101  *              |R2/S2 - f(x)|  < 2**(-61.52)
00102  *
00103  *      5. For inf > x >= 28
00104  *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
00105  *              erfc(x) = tiny*tiny (raise underflow) if x > 0
00106  *                      = 2 - tiny if x<0
00107  *
00108  *      7. Special case:
00109  *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
00110  *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, 
00111  *              erfc/erf(NaN) is NaN
00112  */
00113 
00114 #include <math.h>
00115 
00116 extern "C" {
00117 
00118 /*  assume 32 bit int  */
00119 
00120 typedef int int32_t;
00121 typedef unsigned int u_int32_t;
00122 
00123 /*  assume little endian  */
00124 typedef union 
00125 {
00126   double value;
00127   struct 
00128   {
00129     u_int32_t lsw;
00130     u_int32_t msw;
00131   } parts;
00132 } ieee_double_shape_type;
00133 
00134 
00135 /* Get the more significant 32 bit int from a double.  */
00136 
00137 #define GET_HIGH_WORD(i,d)                                      \
00138 do {                                                            \
00139   ieee_double_shape_type gh_u;                                  \
00140   gh_u.value = (d);                                             \
00141   (i) = gh_u.parts.msw;                                         \
00142 } while (0)
00143 
00144 
00145 /* Set the less significant 32 bits of a double from an int.  */
00146 
00147 #define SET_LOW_WORD(d,v)                                       \
00148 do {                                                            \
00149   ieee_double_shape_type sl_u;                                  \
00150   sl_u.value = (d);                                             \
00151   sl_u.parts.lsw = (v);                                         \
00152   (d) = sl_u.value;                                             \
00153 } while (0)
00154 
00155 
00156 /* Eliminate reference to internal OpenBSD call  */
00157 
00158 #define __ieee754_exp(X) exp(X)
00159 
00160 
00161 static const double
00162 tiny        = 1e-300,
00163 half=  5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
00164 one =  1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
00165 two =  2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
00166         /* c = (float)0.84506291151 */
00167 erx =  8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
00168 /*
00169  * Coefficients for approximation to  erf on [0,0.84375]
00170  */
00171 efx =  1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
00172 efx8=  1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
00173 pp0  =  1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
00174 pp1  = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
00175 pp2  = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
00176 pp3  = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
00177 pp4  = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
00178 qq1  =  3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
00179 qq2  =  6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
00180 qq3  =  5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
00181 qq4  =  1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
00182 qq5  = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
00183 /*
00184  * Coefficients for approximation to  erf  in [0.84375,1.25] 
00185  */
00186 pa0  = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
00187 pa1  =  4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
00188 pa2  = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
00189 pa3  =  3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
00190 pa4  = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
00191 pa5  =  3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
00192 pa6  = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
00193 qa1  =  1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
00194 qa2  =  5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
00195 qa3  =  7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
00196 qa4  =  1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
00197 qa5  =  1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
00198 qa6  =  1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
00199 /*
00200  * Coefficients for approximation to  erfc in [1.25,1/0.35]
00201  */
00202 ra0  = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
00203 ra1  = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
00204 ra2  = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
00205 ra3  = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
00206 ra4  = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
00207 ra5  = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
00208 ra6  = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
00209 ra7  = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
00210 sa1  =  1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
00211 sa2  =  1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
00212 sa3  =  4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
00213 sa4  =  6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
00214 sa5  =  4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
00215 sa6  =  1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
00216 sa7  =  6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
00217 sa8  = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
00218 /*
00219  * Coefficients for approximation to  erfc in [1/.35,28]
00220  */
00221 rb0  = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
00222 rb1  = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
00223 rb2  = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
00224 rb3  = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
00225 rb4  = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
00226 rb5  = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
00227 rb6  = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
00228 sb1  =  3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
00229 sb2  =  3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
00230 sb3  =  1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
00231 sb4  =  3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
00232 sb5  =  2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
00233 sb6  =  4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
00234 sb7  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
00235 
00236         double erf(double x) 
00237 {
00238         int32_t hx,ix,i;
00239         double R,S,P,Q,s,y,z,r;
00240         GET_HIGH_WORD(hx,x);
00241         ix = hx&0x7fffffff;
00242         if(ix>=0x7ff00000) {            /* erf(nan)=nan */
00243             i = ((u_int32_t)hx>>31)<<1;
00244             return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
00245         }
00246 
00247         if(ix < 0x3feb0000) {           /* |x|<0.84375 */
00248             if(ix < 0x3e300000) {       /* |x|<2**-28 */
00249                 if (ix < 0x00800000) 
00250                     return 0.125*(8.0*x+efx8*x);  /*avoid underflow */
00251                 return x + efx*x;
00252             }
00253             z = x*x;
00254             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
00255             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
00256             y = r/s;
00257             return x + x*y;
00258         }
00259         if(ix < 0x3ff40000) {           /* 0.84375 <= |x| < 1.25 */
00260             s = fabs(x)-one;
00261             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
00262             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
00263             if(hx>=0) return erx + P/Q; else return -erx - P/Q;
00264         }
00265         if (ix >= 0x40180000) {         /* inf>|x|>=6 */
00266             if(hx>=0) return one-tiny; else return tiny-one;
00267         }
00268         x = fabs(x);
00269         s = one/(x*x);
00270         if(ix< 0x4006DB6E) {    /* |x| < 1/0.35 */
00271             R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
00272                                 ra5+s*(ra6+s*ra7))))));
00273             S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
00274                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
00275         } else {        /* |x| >= 1/0.35 */
00276             R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
00277                                 rb5+s*rb6)))));
00278             S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
00279                                 sb5+s*(sb6+s*sb7))))));
00280         }
00281         z  = x;  
00282         SET_LOW_WORD(z,0);
00283         r  =  __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
00284         if(hx>=0) return one-r/x; else return  r/x-one;
00285 }
00286 
00287         double erfc(double x) 
00288 {
00289         int32_t hx,ix;
00290         double R,S,P,Q,s,y,z,r;
00291         GET_HIGH_WORD(hx,x);
00292         ix = hx&0x7fffffff;
00293         if(ix>=0x7ff00000) {                    /* erfc(nan)=nan */
00294                                                 /* erfc(+-inf)=0,2 */
00295             return (double)(((u_int32_t)hx>>31)<<1)+one/x;
00296         }
00297 
00298         if(ix < 0x3feb0000) {           /* |x|<0.84375 */
00299             if(ix < 0x3c700000)         /* |x|<2**-56 */
00300                 return one-x;
00301             z = x*x;
00302             r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
00303             s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
00304             y = r/s;
00305             if(hx < 0x3fd00000) {       /* x<1/4 */
00306                 return one-(x+x*y);
00307             } else {
00308                 r = x*y;
00309                 r += (x-half);
00310                 return half - r ;
00311             }
00312         }
00313         if(ix < 0x3ff40000) {           /* 0.84375 <= |x| < 1.25 */
00314             s = fabs(x)-one;
00315             P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
00316             Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
00317             if(hx>=0) {
00318                 z  = one-erx; return z - P/Q; 
00319             } else {
00320                 z = erx+P/Q; return one+z;
00321             }
00322         }
00323         if (ix < 0x403c0000) {          /* |x|<28 */
00324             x = fabs(x);
00325             s = one/(x*x);
00326             if(ix< 0x4006DB6D) {        /* |x| < 1/.35 ~ 2.857143*/
00327                 R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
00328                                 ra5+s*(ra6+s*ra7))))));
00329                 S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
00330                                 sa5+s*(sa6+s*(sa7+s*sa8)))))));
00331             } else {                    /* |x| >= 1/.35 ~ 2.857143 */
00332                 if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
00333                 R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
00334                                 rb5+s*rb6)))));
00335                 S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
00336                                 sb5+s*(sb6+s*sb7))))));
00337             }
00338             z  = x;
00339             SET_LOW_WORD(z,0);
00340             r  =  __ieee754_exp(-z*z-0.5625)*
00341                         __ieee754_exp((z-x)*(z+x)+R/S);
00342             if(hx>0) return r/x; else return two-r/x;
00343         } else {
00344             if(hx>0) return tiny*tiny; else return two-tiny;
00345         }
00346 }
00347 
00348 }
00349 
00350 #else  /* WIN32 */
00351 
00352 int dummy_erf(int i) { return i; }  /* avoid empty translation unit */
00353 
00354 #endif  /* WIN32 */
00355

No veo la necesidad del código, está ya programado en la librería estandar.

Si programas la función error para adaptarla y conseguir diferentes distribuciones de probabilidad, vendría bien ese el aporte al foro.

Un saludo.

PicMinor · « **Respuesta #7 en:** 27 de Mayo de 2014, 04:15:40 »

¡ Saludos al Foro !

Para Picuino:

Dices que no ves la necesidad del código ya que está incluido en la librería standard. Yo estoy trabajando en Visual Basic 6 y no lleva nada incluido.
La necesidad de exactitud con tantos decimales viene dada por unas discrepancias que he encontrado en la acumulación de pequeñas probabilidades respecto a usar un tipo de cálculo u otro.

Finalmente he encontrado dos métodos bastante exactos:

Método de Bevington

Código: [Seleccionar]


Const DOS_DIV_RAIZ_PI = 1.12837916709551
Const RAIZ_2_DIV_2 = 0.707106781186548

' --------------------------------------------------------------------------------------------
' Método de Bevington
' --------------------------------------------------------------------------------------------
Private Function Probabilidad(X As Double, Media As Double, Sigma As Double) As Double
    Dim Z As Double
    Dim Integral As Double
    
    Z = Abs(X - Media) / Sigma
    Integral = IntGaussBevington(Z)                         ' Integral es la Integral entre -Z y +Z
   
    If X > Media Then
        Probabilidad = 0.5 + (Integral / 2)
    Else
        Probabilidad = 0.5 - (Integral / 2)
    End If

End Function

' --------------------------------------------------------------------------------------------
' Método de Bevington
' Devuelve la integral entre -Z y +Z
' --------------------------------------------------------------------------------------------
Private Function IntGaussBevington(Z As Double) As Double
    Dim Termino As Double
    Dim Suma As Double
    Dim Y2 As Double
    Dim Denom As Double
    
    IntGaussBevington = 0
    If Z = 0 Then Exit Function

    Termino = RAIZ_2_DIV_2 * Z
    
    Suma = Termino
    Denom = 1
    Y2 = Z * Z
        
    Do
        Denom = Denom + 2
        Termino = Termino * Y2 / Denom
        Suma = Suma + Termino
        If Termino > 1E+300 Then Exit Do
    Loop Until ((Termino / Suma)) <= 1E-20
        
    IntGaussBevington = DOS_DIV_RAIZ_PI * Suma * Exp(-Y2 / 2)

End Function

Este método tiene el inconveniente de que a veces desborda (Cuando Z es muy alto), de ahí la necesidad de la línea If Termino > 1E+300 Then Exit Do

Otro método que he localizado bastante exacto es el método CumNorm:

Código: [Seleccionar]


' --------------------------------------------------------------------------------------------
' Método CumNorm
' --------------------------------------------------------------------------------------------
Private Function CumNorm(X As Double, Media As Double, Sigma As Double) As Double
    Dim Z As Double
    Dim Exponential As Double
    Dim PolNum As Double
    Dim PolDenom As Double
    
    Z = Abs(X - Media) / Sigma
    
    If Z > 37 Then
        CumNorm = 0
    Else
        Exponential = Exp(-Z * Z / 2)
        If Z < 7.07106781186547 Then
            PolNum = 3.52624965998911E-02 * Z + 0.700383064443688
            PolNum = PolNum * Z + 6.37396220353165
            PolNum = PolNum * Z + 33.912866078383
            PolNum = PolNum * Z + 112.079291497871
            PolNum = PolNum * Z + 221.213596169931
            PolNum = PolNum * Z + 220.206867912376
            PolNum = Exponential * PolNum                             ' Polinomio del Numerador
            
            PolDenom = 8.83883476483184E-02 * Z + 1.75566716318264
            PolDenom = PolDenom * Z + 16.064177579207
            PolDenom = PolDenom * Z + 86.7807322029461
            PolDenom = PolDenom * Z + 296.564248779674
            PolDenom = PolDenom * Z + 637.333633378831
            PolDenom = PolDenom * Z + 793.826512519948
            PolDenom = PolDenom * Z + 440.413735824752                ' Polinomio del Denominador
            
            CumNorm = PolNum / PolDenom
        Else
            PolDenom = Z + 0.65
            PolDenom = Z + 4 / PolDenom
            PolDenom = Z + 3 / PolDenom
            PolDenom = Z + 2 / PolDenom
            PolDenom = Z + 1 / PolDenom
            CumNorm = Exponential / PolDenom / 2.506628274631
        End If
    End If
    
    If (X - Media) > 0 Then CumNorm = 1 - CumNorm
End Function

Y finalmente, por si a alguien le interesa, esta es la implementación de la aproximación de Hastings, que es bastante menos exacto que los dos anteriores:

Código: [Seleccionar]

Const UNO_DIV_RAIZ_2PI = 0.398942280401433

' --------------------------------------------------------------------------------------------
' Método de Hastings
' --------------------------------------------------------------------------------------------
Private Function Hastings(X As Double, Media As Double, Sigma As Double) As Double
    Dim B1 As Double, B2 As Double, B3 As Double, B4 As Double, B5 As Double
    Dim P As Double, S As Double, Z As Double
    
    Z = (X - Media) / Sigma
    
    If Z < -15 Then
        S = Ndf(Z) / Sqr(1 + Z * Z)
    ElseIf Z > 15 Then
        S = 1 - Hastings(-Z, Media, Sigma)
    Else
        P = 0.2316419
        B1 = 0.31938153
        B2 = -0.356563782
        B3 = 1.781477937
        B4 = -1.821255978
        B5 = 1.330274429
        
        S = 1 / (1 + P * Abs(Z))
        S = 1 - Ndf(Z) * (B1 * S + _
                          B2 * S * S + _
                          B3 * S * S * S + _
                          B4 * S * S * S * S + _
                          B5 * S * S * S * S * S)
        
        If (Z <= 0) Then
            S = 1 - S
        End If
        Hastings = S
    End If
End Function

Private Function Ndf(X As Double) As Double
    Ndf = UNO_DIV_RAIZ_2PI * Exp(-X * X * 0.5)              ' 0.398942280401433 = 1/SQR(2*Pi)
End Function

Picuino · « **Respuesta #8 en:** 27 de Mayo de 2014, 08:23:12 »

Tienes razón en que Visual Basic no tiene definida la función de error.

La relación entre la función error y la distribución normal es esta:
(la escribo primero en python por claridad)

Código: Python

def Normal(X, Media, Sigma):
    Z = abs(X - Media) / (math.sqrt(2.0) * Sigma)
    if (X > Media):
       return 0.5 + math.erf(Z)/2.0
    else:
       return 0.5 - math.erf(Z)/2.0

Código: Visual Basic

Function Normal(X As Double, Media As Double, Sigma As Double) As Double
    Dim Z As Double
 
    Z = Abs(X - Media) / Sigma
     
    If X > Media Then
        Normal = 0.5 + Erf(Z) / 2
    Else
        Normal = 0.5 - Erf(Z) / 2
    End If
    
End Function

En el segundo caso, debe definirse en Visual Basic la función de error. Esto es lo que hace el algoritmo de Bevington con la función IntGaussBevington()

Saludos.

Picuino · « **Respuesta #9 en:** 27 de Mayo de 2014, 08:36:08 »

Comparando entre los diferentes algoritmos:

Bevington:
Siempre da valores menores que la función erf de la librería de python.
La precisión relativa es de 15 decimales en el peor de los casos:

Error max = 0
Error min = -1.55431223448e-15
Desviación estandar del error relativo = 8.43866412869e-16

La carga computacional se hace muy grande (hasta 70 iteraciones) en el caso de valores altos de Z

CumNormal:
Tiene error que oscila arriba y abajo de la función erf de la librería de python.
La precisión relativa es de 16 decimales en el peor de los casos:

Error max = 3.83324787322e-16
Error min = -1.60561575907e-16
Desviación estandar del error relativo = 7.18281787181e-17

Está basado en el cálculo
La carga computacional es constante y se reduce a la mitad para valores de Z mayores de 7.07

Este parece ser el mejor método

Hastings:
Este método tiene una precisión menor (7 decimales) y la ventaja de ser más rápido (necesita menos cálculos).

Error max = 1.09166844999e-07
Error min = -1.03855242036e-07
Desviación estandar del error relativo = 3.98861416192e-08

Los cálculos están realizados en Python 2.7 con Z en el intervalo de 0 a 7 en pasos de 0.1

El código es:

Código: Python

# -*- coding: cp1252 -*-
#
# Comparativa entre métodos de calcular la función de distribución Normal
#
 
import math
 
def main():
    print "\n\nCompare CumNorm with erf:"
    compare(CumNorm)
 
    print "\n\nCompare Bevington with erf:"
    compare(Bevington)
 
    print "\n\nCompare Hastings with erf:"
    compare(Hastings)
 
    #raw_input("Pulsa Enter")
 
 
def compare(func):
    xr = [x/10.0 for x in range(0, 8*10)]
    maxerr = minerr = desvest = n = 0
    for x in xr:
        y1 = lib_erf(x, 0, 1)
        y2 = func(x, 0, 1)
        err = error_rel(y2, y1)
        desvest = desvest + err**2
        n = n + 1
        if (err > maxerr): maxerr = err
        if (err < minerr): minerr = err
        print "  x=%1.1f\ty=%1.16f\tError=%-1.3e" % (x, y2, error_rel(y2, y1))
    print "Error max =", maxerr
    print "Error min =", minerr
    print "Desviación estandar =", math.sqrt(desvest/n)
 
def error_rel(x, ref):
    if (ref == 0):
        return 0
    return (x-ref)/ref
 
 
# --------------------------------------------------------------------------------------------
#  Método de la función de error de la librería estandar
# --------------------------------------------------------------------------------------------
def lib_erf(X, Media, Sigma):
    Z = abs(X - Media) / (math.sqrt(2.0) * Sigma)
    if (X > Media):
       return 0.5 + math.erf(Z)/2.0
    else:
       return 0.5 - math.erf(Z)/2.0
 
 
# --------------------------------------------------------------------------------------------
#  Método de la función de error CumNorm
# --------------------------------------------------------------------------------------------
def CumNorm(X, Media, Sigma):
 
    Z = abs(X - Media) / Sigma
 
    if (Z > 37):
        return 0
    
    else:
        Exponential = math.exp(-Z * Z / 2)
        
        if (Z < 7.07106781186547):
            PolNum = 3.52624965998911E-02 * Z + 0.700383064443688
            PolNum = PolNum * Z + 6.37396220353165
            PolNum = PolNum * Z + 33.912866078383
            PolNum = PolNum * Z + 112.079291497871
            PolNum = PolNum * Z + 221.213596169931
            PolNum = PolNum * Z + 220.206867912376
            PolNum = Exponential * PolNum
            
            PolDenom = 8.83883476483184E-02 * Z + 1.75566716318264
            PolDenom = PolDenom * Z + 16.064177579207
            PolDenom = PolDenom * Z + 86.7807322029461
            PolDenom = PolDenom * Z + 296.564248779674
            PolDenom = PolDenom * Z + 637.333633378831
            PolDenom = PolDenom * Z + 793.826512519948
            PolDenom = PolDenom * Z + 440.413735824752
            
            valor = PolNum / PolDenom
 
        else:
            PolDenom = Z + 0.65
            PolDenom = Z + 4 / PolDenom
            PolDenom = Z + 3 / PolDenom
            PolDenom = Z + 2 / PolDenom
            PolDenom = Z + 1 / PolDenom
            valor = Exponential / PolDenom / 2.506628274631
    
    if (X > 0):
        valor = 1 - valor
 
    return valor
 
 
# --------------------------------------------------------------------------------------------
#  Método de la función error de Bevington
# --------------------------------------------------------------------------------------------
def Bevington(X, Media, Sigma):
 
    DOS_DIV_RAIZ_PI = 1.12837916709551
    RAIZ_2_DIV_2 = 0.707106781186548
 
    def erf_Bevington(Z):
        if (Z == 0):
           return 0
        Termino = RAIZ_2_DIV_2 * Z
        Suma = Termino
        Denom = 1
        Y2 = Z * Z
        N = 0
        while(True):
            N = N + 1
            Denom = Denom + 2
            Termino = Termino * Y2 / Denom
            Suma = Suma + Termino
            if (Termino > 1E+300):
                break
            if ((Termino / Suma) <= 1E-20):
                break
        return DOS_DIV_RAIZ_PI * Suma * math.exp(-Y2 / 2)
 
    Z = abs(X - Media) / Sigma
   
    if (X > Media):
        return 0.5 + (erf_Bevington(Z) / 2)
    else:
        return 0.5 - (erf_Bevington(Z) / 2)
 
 
 
 
# --------------------------------------------------------------------------------------------
#  Método de Hastings
# --------------------------------------------------------------------------------------------
def Hastings(X, Media, Sigma):
    
    def Ndf(X):
       UNO_DIV_RAIZ_2PI = 0.398942280401433
       return UNO_DIV_RAIZ_2PI * math.exp(-X * X * 0.5)
    
    Z = (X - Media) / Sigma
    
    if (Z < -15):
        S = Ndf(Z) / math.sqrt(1 + Z * Z)
    elif (Z > 15):
        S = 1 - Hastings(-Z, Media, Sigma)
    else:
        P = 0.2316419
        B1 = 0.31938153
        B2 = -0.356563782
        B3 = 1.781477937
        B4 = -1.821255978
        B5 = 1.330274429
        
        S = 1 / (1 + P * abs(Z))
        S = 1 - Ndf(Z) * S*(B1 + S*(B2 + S*(B3 + S*(B4 +  S*B5))))
        
        if (Z <= 0):
            S = 1 - S
 
    return S

Saludos.

Picuino · « **Respuesta #10 en:** 27 de Mayo de 2014, 18:09:29 »

Otro método, en C:

Código: C

/*
   C program for floating point error function
 
   erf(x) returns the error function of its argument
   erfc(x) returns 1.0-erf(x)
 
   erf(x) is defined by
   ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
 
   the entry for erfc is provided because of the
   extreme loss of relative accuracy if erf(x) is
   called for large x and the result subtracted
   from 1. (e.g. for x= 10, 12 places are lost).
 
   There are no error returns.
   Coefficients for large x are #5667 from Hart & Cheney (18.72D).
*/
#include <errno.h>
#include "cmath.h"
 
#define M 7
#define N 9
 
static const double torp = 1.1283791670955125738961589031;
static const double p1[] = {
    0.804373630960840172832162e5,
    0.740407142710151470082064e4,
    0.301782788536507577809226e4,
    0.380140318123903008244444e2,
    0.143383842191748205576712e2,
   -0.288805137207594084924010e0,
    0.007547728033418631287834e0,
};
static const double q1[]  = {
    0.804373630960840172826266e5,
    0.342165257924628539769006e5,
    0.637960017324428279487120e4,
    0.658070155459240506326937e3,
    0.380190713951939403753468e2,
    0.100000000000000000000000e1,
    0.0,
};
static const double p2[]  = {
    0.1826334884229511259216899900e4,
    0.2898029329216765561127584600e4,
    0.2320439590251635247384768711e4,
    0.1143262070703886173606073338e4,
    0.3685196154710010637133875746e3,
    0.7708161730368428609781633646e2,
    0.9675807882987265400604202961e1,
    0.5641877825507397413087057563e0,
    0.0,
};
static const double q2[]  = {
    0.1826334884229511259557643800e4,
    0.4958827564721140714954384220e4,
    0.6089542423272443550463306800e4,
    0.4429612803883682726711528526e4,
    0.2094384367789539593790281779e4,
    0.6617361207107653469211984771e3,
    0.1371255960500622202878443578e3,
    0.1714980943627607849376131193e2,
    1.0,
};
 
double erf(double arg)
{
   int sign, i;
   double argsq, d, n;
 
   sign = 1;
   if(arg < 0.0)
   {
      arg = -arg;
      sign = -1;
   }
   if(arg < 0.5)
   {
      argsq = arg*arg;
      for(n=0,d=0,i=M-1; i>=0; i--)
      {
         n = n*argsq + p1[i];
         d = d*argsq + q1[i];
      }
      if (sign<0) return(-torp*arg*n/d);
      return(torp*arg*n/d);
   }
   if(arg >= 10.0)  return(sign*1.);
   if (sign<0) return(-1.0 + erfc(arg));
   return(1.0 - erfc(arg));
}
 
double erfc(double arg)
{
   double n, d;
   int i;
 
   errno = 0;
   if(arg < 0.0) return(2.0 - erfc(-arg));
   if(arg < 0.5) return(1.0 - erf(arg));
   if(arg >= 10.0) return(0.0);
 
   for(n=0,d=0,i=N-1; i>=0; i--)
   {
      n = n*arg + p2[i];
      d = d*arg + q2[i];
   }
   return(exp(-arg*arg)*n/d);
}

Las rutinas tambien incluyen la función erfc que es:

erfc = 1 - erf()

Conviene calcular erfc por separado, porque la fórmula anterior da un error muy grande para valores altos de Z.

Un saludo.

Picuino · « **Respuesta #11 en:** 27 de Mayo de 2014, 18:19:32 »

Los coeficientes vienen del libro:

Computer Approximations. John F. Hart & Cheney. Ed. John Wiley & Sons

Un saludo.

PicMinor · « **Respuesta #12 en:** 28 de Mayo de 2014, 03:12:16 »

¡ Saludos al foro!

Para Picuino:

Gracias por tu excelente estudio de los algoritmos. Para el trabajo que estoy haciendo he decidido usar el método cumNorm. Si me queda tiempo intentaré traducir la rutina en C a Visual Basic. ¿Esta rutina es todavía más precisa que el método cumNorm?

Picuino · « **Respuesta #13 en:** 28 de Mayo de 2014, 04:24:41 »

No lo creo. CumNorm llega a 16 decimales de precisión, que es el límite de los números en coma flotante double (de 64 bits, con 54 bits de mantisa).

Donde si que puedes ganar precisión es en la función erfc para valores altos de X. Para eso está diseñada.

Un saludo.

Picuino · « **Respuesta #14 en:** 28 de Mayo de 2014, 06:21:49 »

Comparativa de varios métodos para calcular la Distribución Normal
La mayoría están basados en calcular la función de error erf()

Código para la comparación de métodos:

Código: Python

# -*- coding: cp1252 -*-
 
import math
 
def main():
    print "\n\nCompare CumNorm   with system erf:"
    compare(CumNorm)
 
    print "\n\nCompare Bevington with system erf:"
    compare(Bevington)
 
    print "\n\nCompare Hastings  with system erf:"
    compare(Hastings)
 
    print "\n\nCompare Hart      with system erf:"
    compare(Norm_Hart)
 
    #raw_input("Pulsa Enter")
 
 
def compare(func):
    xr = [x/10.0 for x in range(0, 8*10)]
    maxerr = minerr = desvest = n = 0
    for x in xr:
        y1 = normalize_erf(math.erf, x, 0, 1)
        y2 = func(x, 0, 1)
        err = error_rel(y2, y1)
        desvest = desvest + err**2
        n = n + 1
        if (err > maxerr): maxerr = err
        if (err < minerr): minerr = err
        # print "  x=%1.1f\ty=%1.16f\tError=%-1.3e" % (x, y2, error_rel(y2, y1))
    print "   Error max =", maxerr
    print "   Error min =", minerr
    print "   Desviación estandar =", math.sqrt(desvest/n)
 
 
def error_rel(x, ref):
    if (ref == 0):
        return 0
    return (x-ref)/ref
 
 
 
# --------------------------------------------------------------------------------------------
#  Normalizacion de la funcion de error
# --------------------------------------------------------------------------------------------
def normalize_erf(erf, X, Media, Sigma):
    Z = abs(X - Media) / (math.sqrt(2.0) * Sigma)
    if (X > Media):
       return 0.5 + erf(Z)/2.0
    else:
       return 0.5 - erf(Z)/2.0
 
 
# --------------------------------------------------------------------------------------------
#  Método de Hart and Cheney
# --------------------------------------------------------------------------------------------    
def Norm_Hart(X, Media, Sigma):
    return normalize_erf(erf_Hart, X, Media, Sigma)
 
 
def erf_Hart(arg):
   torp = 1.1283791670955125738961589031 
   sign = 1
   if(arg < 0.0):
      arg = -arg
      sign = -1
     
   if(arg < 0.5):
      argsq = arg*arg
      n = 0.007547728033418631287834e0
      n = n*argsq - 0.288805137207594084924010e0
      n = n*argsq + 0.143383842191748205576712e2
      n = n*argsq + 0.380140318123903008244444e2
      n = n*argsq + 0.301782788536507577809226e4
      n = n*argsq + 0.740407142710151470082064e4
      n = n*argsq + 0.804373630960840172832162e5
      d = 0.100000000000000000000000e1
      d = d*argsq + 0.380190713951939403753468e2
      d = d*argsq + 0.658070155459240506326937e3
      d = d*argsq + 0.637960017324428279487120e4
      d = d*argsq + 0.342165257924628539769006e5
      d = d*argsq + 0.804373630960840172826266e5
      
      if (sign<0):
         return (-torp*arg*n/d)
      else:
         return (torp*arg*n/d)
 
   if (arg >= 10.0):
      return sign*1.0
    
   if (sign < 0):
      return (-1.0 + erfc_Hart(arg))
   else:
      return(1.0 - erfc_Hart(arg));
 
 
 
def erfc_Hart(arg):
   errno = 0
   if (arg < 0.0):
      return(2.0 - erfc(-arg))
   if (arg < 0.5):
      return(1.0 - erf(arg))
   if (arg >= 10.0):
      return(0.0)
   
   n = 0.5641877825507397413087057563e0
   n = n*arg + 0.9675807882987265400604202961e1
   n = n*arg + 0.7708161730368428609781633646e2
   n = n*arg + 0.3685196154710010637133875746e3
   n = n*arg + 0.1143262070703886173606073338e4
   n = n*arg + 0.2320439590251635247384768711e4
   n = n*arg + 0.2898029329216765561127584600e4
   n = n*arg + 0.1826334884229511259216899900e4
   d = 1.0
   d = d*arg + 0.1714980943627607849376131193e2
   d = d*arg + 0.1371255960500622202878443578e3
   d = d*arg + 0.6617361207107653469211984771e3
   d = d*arg + 0.2094384367789539593790281779e4
   d = d*arg + 0.4429612803883682726711528526e4
   d = d*arg + 0.6089542423272443550463306800e4
   d = d*arg + 0.4958827564721140714954384220e4
   d = d*arg + 0.1826334884229511259557643800e4
   return (math.exp(-arg*arg)*n/d);
 
 
 
# --------------------------------------------------------------------------------------------
#  Método de la función de error CumNorm
# --------------------------------------------------------------------------------------------
def CumNorm(X, Media, Sigma):
 
    Z = abs(X - Media) / Sigma
 
    if (Z > 37):
        return 0
    
    else:
        Exponential = math.exp(-Z * Z / 2)
        
        if (Z < 7.07106781186547):
            PolNum = 3.52624965998911E-02 * Z + 0.700383064443688
            PolNum = PolNum * Z + 6.37396220353165
            PolNum = PolNum * Z + 33.912866078383
            PolNum = PolNum * Z + 112.079291497871
            PolNum = PolNum * Z + 221.213596169931
            PolNum = PolNum * Z + 220.206867912376
            PolNum = Exponential * PolNum
            
            PolDenom = 8.83883476483184E-02 * Z + 1.75566716318264
            PolDenom = PolDenom * Z + 16.064177579207
            PolDenom = PolDenom * Z + 86.7807322029461
            PolDenom = PolDenom * Z + 296.564248779674
            PolDenom = PolDenom * Z + 637.333633378831
            PolDenom = PolDenom * Z + 793.826512519948
            PolDenom = PolDenom * Z + 440.413735824752
            
            valor = PolNum / PolDenom
 
        else:
            PolDenom = Z + 0.65
            PolDenom = Z + 4 / PolDenom
            PolDenom = Z + 3 / PolDenom
            PolDenom = Z + 2 / PolDenom
            PolDenom = Z + 1 / PolDenom
            valor = Exponential / PolDenom / 2.506628274631
    
    if (X > 0):
        valor = 1 - valor
 
    return valor
 
 
# --------------------------------------------------------------------------------------------
#  Método de la función error de Bevington
# --------------------------------------------------------------------------------------------
def Bevington(X, Media, Sigma):
 
    DOS_DIV_RAIZ_PI = 1.12837916709551
    RAIZ_2_DIV_2 = 0.707106781186548
 
    def erf_Bevington(Z):
        if (Z == 0):
           return 0
        Termino = RAIZ_2_DIV_2 * Z
        Suma = Termino
        Denom = 1
        Y2 = Z * Z
        N = 0
        while(True):
            N = N + 1
            Denom = Denom + 2
            Termino = Termino * Y2 / Denom
            Suma = Suma + Termino
            if (Termino > 1E+300):
                break
            if ((Termino / Suma) <= 1E-20):
                break
        return DOS_DIV_RAIZ_PI * Suma * math.exp(-Y2 / 2)
 
    Z = abs(X - Media) / Sigma
   
    if (X > Media):
        return 0.5 + (erf_Bevington(Z) / 2)
    else:
        return 0.5 - (erf_Bevington(Z) / 2)
 
 
 
 
# --------------------------------------------------------------------------------------------
#  Método de Hastings
# --------------------------------------------------------------------------------------------
def Hastings(X, Media, Sigma):
    
    def Ndf(X):
       UNO_DIV_RAIZ_2PI = 0.398942280401433
       return UNO_DIV_RAIZ_2PI * math.exp(-X * X * 0.5)
    
    Z = (X - Media) / Sigma
    
    if (Z < -15):
        S = Ndf(Z) / math.sqrt(1 + Z * Z)
    elif (Z > 15):
        S = 1 - Hastings(-Z, Media, Sigma)
    else:
        P = 0.2316419
        B1 = 0.31938153
        B2 = -0.356563782
        B3 = 1.781477937
        B4 = -1.821255978
        B5 = 1.330274429
        
        S = 1 / (1 + P * abs(Z))
        S = 1 - Ndf(Z) * S*(B1 + S*(B2 + S*(B3 + S*(B4 +  S*B5))))
        
        if (Z <= 0):
            S = 1 - S
 
    return S
 
 
main()

Resultados:
Compare CumNorm with system erf: Error max = 3.83324787322e-16 Error min = -1.60561575907e-16 Desviación estandar = 7.18281787181e-17 Compare Bevington with system erf: Error max = 0 Error min = -1.55431223448e-15 Desviación estandar = 8.43866412869e-16 Compare Hastings with system erf: Error max = 1.09166844999e-07 Error min = -1.03855242036e-07 Desviación estandar = 3.98861416192e-08 Compare Hart with system erf: Error max = 1.22921127135e-16 Error min = -1.31958157439e-16 Desviación estandar = 3.97910289557e-17

El método basado en los coeficientes de Hart y Cheney es el que más se parece a la función de error del sistema aportada por Python

Saludos.

TODOPIC

Autor Tema: Integración de la Campana de Gauss (Leído 8721 veces)

PicMinor

Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

PicMinor

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

PicMinor

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

PicMinor

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss

Picuino

Re: Integración de la Campana de Gauss