Интеграл от функции нормального распределения

[Puma]

Весь интернет обыскал, нигде не нашел. Решил сам написать. Вот выкладываю думаю кому-нибудь пригодиться!


/*************************************************
Вычисление функции erf

Абсолютная погрешность полученного значения не превосходит 10^-8.

Входные параметры:
    x   -   аргумент

*************************************************/

function erf($x)
{

    if( $x<0.5 )
    {
        $xsq = $x*$x;
        $p = 0.007547728033418631287834;
        $p = 0.288805137207594084924010+$xsq*$p;
        $p = 14.3383842191748205576712+$xsq*$p;
        $p = 38.0140318123903008244444+$xsq*$p;
        $p = 3017.82788536507577809226+$xsq*$p;
        $p = 7404.07142710151470082064+$xsq*$p;
        $p = 80437.3630960840172832162+$xsq*$p;
        $q = 0.0;
        $q = 1.00000000000000000000000+$xsq*$q;
        $q = 38.0190713951939403753468+$xsq*$q;
        $q = 658.070155459240506326937+$xsq*$q;
        $q = 6379.60017324428279487120+$xsq*$q;
        $q = 34216.5257924628539769006+$xsq*$q;
        $q = 80437.3630960840172826266+$xsq*$q;
        $result = s*1.1283791670955125738961589031*$x*$p/$q;
        return $result;
    }

    if( $x>=10 )
    {
        $result = $x;
        return $result;
    }

    $result = $x*(1-erfc($x));
    return $result;
}


/*************************************************Вычисление дополнительной функции erfc

Абсолютная погрешность полученного значения не превосходит 10^-8.

Входные параметры:
    x   -   аргумент

*************************************************/

function erfc($x)
{
    if( $x<0 )
    {
        $result = 2-erfc(-$x);
        return $result;
    }
    if( $x<0.5 )
    {
        $result = 1.0-erf($x);
        return $result;
    }
    if( $x>=10 )
    {
        $result = 0;
        return $result;
    }
    $p = 0.0;
    $p = 0.5641877825507397413087057563+$x*$p;
    $p = 9.675807882987265400604202961+$x*$p;
    $p = 77.08161730368428609781633646+$x*$p;
    $p = 368.5196154710010637133875746+$x*$p;
    $p = 1143.262070703886173606073338+$x*$p;
    $p = 2320.439590251635247384768711+$x*$p;
    $p = 2898.0293292167655611275846+$x*$p;
    $p = 1826.3348842295112592168999+$x*$p;
    $q = 1.0;
    $q = 17.14980943627607849376131193+$x*$q;
    $q = 137.1255960500622202878443578+$x*$q;
    $q = 661.7361207107653469211984771+$x*$q;
    $q = 2094.384367789539593790281779+$x*$q;
    $q = 4429.612803883682726711528526+$x*$q;
    $q = 6089.5424232724435504633068+$x*$q;
    $q = 4958.82756472114071495438422+$x*$q;
    $q = 1826.3348842295112595576438+$x*$q;
    $result = exp(-sqrt($x))*$p/$q;
    return $result;
}

/*************************************************Вычисление интегральной функции вероятности нормального распределения

Входные параметры:
    x   -   аргумент

*************************************************/

function normaldistribution($x)
{
    $result = 0.5*(erf($x/1.41421356237309504880)+1);
    return $result;
}