Logo Search packages:      
Sourcecode: openturns version File versions  Download package

DistFunc.cxx

Go to the documentation of this file.
//                                               -*- C++ -*-
/**
 *  @file  DistFunc.cxx
 *  @brief OpenTURNS wrapper to a library of special functions
 *
 *  (C) Copyright 2005-2007 EDF-EADS-Phimeca
 *
 *  This library is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU Lesser General Public
 *  License as published by the Free Software Foundation; either
 *  version 2.1 of the License.
 *
 *  This library is distributed in the hope that it will be useful
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *  Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public
 *  License along with this library; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 *
 *  @author: $LastChangedBy: dutka $
 *  @date:   $LastChangedDate: 2008-10-31 11:52:04 +0100 (ven 31 oct 2008) $
 *  Id:      $Id: DistFunc.cxx 995 2008-10-31 10:52:04Z dutka $
 */
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include "DistFunc.hxx"
#include "RandomGenerator.hxx"
#include "Log.hxx"
#include "Nct.hxx"
#include "SquareMatrix.hxx"

#define USE_DCDFLIB
#ifdef USE_DCDFLIB
#include "Dcdflib.hxx"
#endif

#include "SpecFunc.hxx"

namespace OpenTURNS {

  namespace Uncertainty {

    namespace Distribution {

      typedef Base::Stat::RandomGenerator RandomGenerator;
      typedef Base::Type::SquareMatrix    SquareMatrix;
      typedef Base::Common::Log           Log;

      const NumericalScalar DistFunc::Precision;
      const UnsignedLong    DistFunc::MaximumIteration;
      const UnsignedLong    DistFunc::NumberOfBandNormalZigurrat;
      const NumericalScalar DistFunc::NormalZigguratTail;
      const NumericalScalar DistFunc::NormalZigguratAbscissa[DistFunc::NumberOfBandNormalZigurrat + 1] = {0.0, 0.271599048510693754998163050702, 0.361898145612134821317157669025, 0.425392047157312675983415424660, 0.476131717888617028166988513439, 0.519219013293501476350033678035, 0.557137093511580298809061311694, 0.591298984443313319787098187809, 0.622592649542783642315521009315, 0.651616031063230052452116746597, 0.678792877347767760686380971868, 0.704435546364376654764501920124, 0.728781608504968321114578996184, 0.752016422469955992826848889422, 0.774287710877425340184686032696, 0.795715332287609618751456615712, 0.816398043479172731605919440849, 0.836418306930649053279379059005, 0.855845789056994122927551923075, 0.874739957924955634206040955649, 0.893152046999867275098199151897, 0.911126563328297873672484801050, 0.928702462334730770821130245805, 0.945914074632515687029155431559, 0.962791845646485921191452370469, 0.979362932051085821966553188665, 0.995651687354039595085786832964, 1.01168006070531516232965392900, 1.02746792709172237058780318928, 1.04303336277130264040166296783, 1.05839287662899292124745211169, 1.07356160576994269038108408924, 1.08855348188453823734559909469, 1.10338137356245140332933074704, 1.11805720869053844255517390583, 1.13259208026135637042889232625, 1.14699633828758735801507046149, 1.16127967002039429336339539399, 1.17545117027533417287473048372, 1.18951940335453972086012397336, 1.20349245780083403835096451764, 1.21737799501490321160206827336, 1.23118329260039645151126944958, 1.24491528316597898248595200283, 1.25858058920182308169327058845, 1.27218555455602081689118701416, 1.28573627296019479845973032796, 1.29923861399021398514811910954, 1.31269824679504233415845967597, 1.32612066188248490208347357933, 1.33951119121344617033356972221, 1.35287502682506865937549704692, 1.36621723817679783092212387778, 1.37954278839122803847614338568, 1.39285654954287803582443862287, 1.40616331713229672637906532114, 1.41946782386968059570797384556, 1.43277475288114273351645465423, 1.44608875044162355949759042785, 1.45941443833094318522780953546, 1.47275642590347769759399918804, 1.48611932195724762500809972081, 1.49950774648472030334842211111, 1.51292634238526126953827654531, 1.52637978721786115033208986070, 1.53987280507247501439918234146, 1.55341017863902384055596007811, 1.56699676155482666996780634172, 1.58063749111398211301308134969, 1.59433740142604557618725800299, 1.60810163711632290631602733261, 1.62193546766631570500142829940, 1.63584430250042915164619942820, 1.64983370693414137939082821110, 1.66390941910962125167951454089, 1.67807736805749698528091063361, 1.69234369303839820145325283038, 1.70671476433535309787605152682, 1.72119720568852405198043944381, 1.73579791858759592141174442291, 1.75052410866497944355179477982, 1.76538331446556845629758885967, 1.78038343890695608962355720514, 1.79553278378881937566665179110, 1.81084008776290159625982126988, 1.82631456823722069305589567937, 1.84196596776173086857556181334, 1.85780460553003206309142077930, 1.87384143473579289060235282188, 1.89008810664697886815311451817, 1.90655704241032628459974673316, 1.92326151377851307313333422169, 1.94021573417040550238786161927, 1.95743496173982213662503222153, 1.97493561645224006631535036596, 1.99273541356693511100128687843, 2.01085351641383291996867861457, 2.02931071196545767371544926885, 2.04812961346843841880165860053, 2.06733489536060750602383910558, 2.08695356691825680949271767439, 2.10701529263367409420415792810, 2.12755276932450072216751946158, 2.14860217257314763105233505342, 2.17020368849378343938267967922, 2.19240215131716394363476960255, 2.21524781328299046507842208032, 2.23879728143015418864904200351, 2.26311466694490835351573544076, 2.28827300805157643267084011891, 2.31435604894699812266939486476, 2.34146048795978434043531292659, 2.36969885260647551209738645471, 2.39920322494817055544917817450, 2.43013013978464913728162829644, 2.46266713119979698270547205291, 2.49704164517949055368308433652, 2.53353343078795745430738206406, 2.57249218737695205289488067611, 2.61436340874381837128188016042, 2.65972749227020134812347957976, 2.70936127691136762243924804229, 2.76433956507998809402134093896, 2.82621272644033894046564603318, 2.89734153514829733018851271394, 2.98159419140365968437694926385, 3.08601280427975485670709771297, 3.22573136821851522103605577983, 3.44508288805539135654449538289, 3.71537970891694553086615539964};
      const NumericalScalar DistFunc::NormalZigguratRatio[DistFunc::NumberOfBandNormalZigurrat] = {0., 0.750484775353838553484969850043, 0.850740271310955168539797476445, 0.893433541969631939167066002439, 0.917015181837094471851335493153, 0.931941203234046233441298028733, 0.942225689827803086945282581056, 0.949736532992396229158449917333, 0.955459380775071654971678008833, 0.959962976643589437257667346996, 0.963598275031758611303047864447, 0.966593473467949497522283939333, 0.969103315737874152720557392852, 0.971236417555650662146408519772, 0.973071247290683962507094254337, 0.974665898140297597377304564176, 0.976064293086860561017929573224, 0.977300253883645249170475123945, 0.978400244899316179173188425569, 0.979385269129977846895884043302, 0.980272206900904494815390564396, 0.981074779362329241400634997295, 0.981804253938740083291033345358, 0.982469968882383589521668758583, 0.983079728809119850562796554795, 0.983640106766412054356146072398, 0.984156677615944084392929857459, 0.984634200280007552311796735515, 0.985076761458307086707253212313, 0.985487889991653471429663256532, 0.985870648634019490358385815715, 0.986227708271493343086583402375, 0.986561408382272362985077361984, 0.986873806622761920341336293274, 0.987166719753625799389459567589, 0.987441757619090675912393772106, 0.987700351516051174211224554768, 0.987943778003453440265050048019, 0.988173178983434816520165815368, 0.988389578716739479383699175238, 0.988593898303624926534040096469, 0.988786968058724293246032191279, 0.988969538127397455898436508912, 0.989142287626960404173170844597, 0.989305832545044324977737622638, 0.989460732586336934407985727862, 0.989607497125912170683173044980, 0.989746590400581327272272451975, 0.989878436047901646747803160708, 0.990003421084648827869197604239, 0.990121899401909443020286447879, 0.990234194841858232813880454444, 0.990340603911263982506239462637, 0.990441398178427316016322972125, 0.990536826393283914925721684389, 0.990627116364558512218324337623, 0.990712476622928023181406791844, 0.990793097894984079572175093174, 0.990869154409244076145976128873, 0.990940805052437824691268287474, 0.991008194391705548082954937064, 0.991071453576109215637707296064, 0.991130701128922528307152229819, 0.991186043640474608746669498017, 0.991237576369836074033189634373, 0.991285383762317554082871105898, 0.991329539888568903748275039752, 0.991370108809995454035219189473, 0.991407144874221943150024047055, 0.991440692943413945013043999353, 0.991470788557391074579718311726, 0.991497458032617503483410551001, 0.991520718497315401540555107636, 0.991540577862097828628588285611, 0.991557034724640724341601679445, 0.991570078205989141038468359572, 0.991579687715098995526722137843, 0.991585832637127947522049524329, 0.991588471939779518821684067585, 0.991587553690640507148992013935, 0.991583014476894319513703858007, 0.991574778716995497249215923087, 0.991562757851797002828992213523, 0.991546849400162799407295751574, 0.991526935861188886574791544186, 0.991502883441690831444489168144, 0.991474540583463335923580694623, 0.991441736259812196096950028863, 0.991404278004792265983740672627, 0.991361949631191261607889802980, 0.991314508584238996353362687132, 0.991261682866858262508333044405, 0.991203167458441781253777050367, 0.991138620131902842410876661610, 0.991067656552139074390092217899, 0.990989844511792319927328343055, 0.990904697125582941945109300521, 0.990811664760263111407350456154, 0.990710125420284608160345942834, 0.990599373235399951429113359060, 0.990478604599807611123088277090, 0.990346901385037190355276297326, 0.990203210479193398963460161789, 0.990046318677290513694722835952, 0.989874821637971854574784332163, 0.989687085197041986702889368351, 0.989481196737866229216740827939, 0.989254903487686985183950871399, 0.989005533422740516025904697311, 0.988729892746066823168948592841, 0.988424131369219266125861236354, 0.988083564029399232263953634600, 0.987702428858508882825722518646, 0.987273556123533683467427764224, 0.986787905274353495741507623455, 0.986233904410023291038018271182, 0.985596485459788251551070665927, 0.984855636576794042676975880029, 0.983984161793717454104724720040, 0.982944085941803527491004884902, 0.981680632603656288469381808499, 0.980111601026471711171610577703, 0.978107394117398550145114447318, 0.975450319596402797690182574081, 0.971742413337712360174408714492, 0.966163908091604482763906551568, 0.956686237014236193571274296719, 0.936329102386070277297433094023, 0.927249206800360320095099544455};

      /**************************************************************************************************************************/
      /* Normalized Beta distribution, i.e. with a PDF equals to x ^ (p1 - 1) . (1 - x) ^ (p2 - 1) / Beta(p1, p2) . (0 < x < 1) */
      /**************************************************************************************************************************/
      /* CDF */
      NumericalScalar DistFunc::pBeta(const NumericalScalar p1, const NumericalScalar p2, const NumericalScalar x, const Bool tail)
      {
      return SpecFunc::BetaRatioInc(p1, p2, x);
      }
      /* CDF Inverse */
      NumericalScalar DistFunc::qBeta(const NumericalScalar p1, const NumericalScalar p2, const NumericalScalar x, const Bool tail)
      {
      return SpecFunc::BetaRatioIncInv(p1, p2, x);
      }
      /* Random number generation
       We use the algorithm of Cheng (1978), Johnk, Atkinson and Whittaker (1979) 1 & 2 described in:
       Luc Devroye, "Non-Uniform RandomVariate Generation", Springer-Verlag, 1986, available online at:
       http://cg.scs.carleton.ca/~luc/nonuniformrandomvariates.zip
       and with the important errata at:
       http://cg.scs.carleton.ca/~luc/errors.pdf
      */
      NumericalScalar DistFunc::rBeta(const NumericalScalar p1, const NumericalScalar p2)
      {
      // Strategy:
      // If (a = 1 and b = 1), Beta(1,1) = Uniform(0,1)
      // If (a = 1 or b = 1), analytic cases
      // If (a + b <= 1), Johnk
      // If (a + b > 1):
      // If (a < 1 and b < 1), Atkinson and Whittaker 1
      // If (a < 1 and b > 1) or (a > 1 and b < 1), Atkinson and Whittaker 2
      // If (a > 1 and b > 1) Cheng.
      if ((p1 == 1.0) && (p2 == 1.0)) return RandomGenerator::Generate();
      // Special case for p1 = 1 or p2 = 1
      if (p1 == 1.0)
        {
          return 1.0 - pow(RandomGenerator::Generate(), 1.0 / p2);
        }
      if (p2 == 1.0)
        {
          return pow(RandomGenerator::Generate(), 1.0 / p1);
        }
      // Now, the more general cases
      NumericalScalar minp(std::min(p1, p2));
      NumericalScalar maxp(std::max(p1, p2));
      NumericalScalar sum(p1 + p2);
      if (sum <= 1.0) // Johnk
        {
          // Use of logarithms to avoid underflow if minp << 1 (here 1e-3, quite arbitrary). It is the
          // only case where this kind of trick is useful. It should not be much slower than the usual
          // version using pow()
          if (minp < 1e-3)
            {
            for (;;)
              {
                NumericalScalar u(RandomGenerator::Generate());
                NumericalScalar v(RandomGenerator::Generate());
                NumericalScalar logx(log(u) / p1);
                NumericalScalar logy(log(v) / p2);
                NumericalScalar logsum((logx > logy) ? logx + log(1.0 + exp(logy - logx)) : logy + log(1.0 + exp(logx - logy)));
                // Acceptation step
                if (logsum <= 0.0) return exp(logy - logsum);
              }
            }
          // Usual form of the algorithm
          for (;;)
            {
                NumericalScalar u(RandomGenerator::Generate());
                NumericalScalar v(RandomGenerator::Generate());
                NumericalScalar x(pow(u, 1.0 / p1));
                NumericalScalar y(pow(v, 1.0 / p2));
                // Acceptation step
                if (x + y <= 1.0) return x / (x + y);
            }
        } // End Johnk
      // Now, sum > 1 for all the remaining cases
      if (minp > 1.0) // Cheng
        {
          NumericalScalar lambda(sqrt((sum - 2.0) / (2.0 * p1 * p2 - sum)));
          NumericalScalar c(minp + 1.0 / lambda);
          for (;;)
            {
            NumericalScalar u1(RandomGenerator::Generate());
            NumericalScalar u2(RandomGenerator::Generate());
            NumericalScalar v(lambda * log(u1 / (1.0 - u1)));
            NumericalScalar w(minp * exp(v));
            NumericalScalar z(u1 * u1 * u2);
            // 1.386294361119890618834464 = log(4)
            NumericalScalar r(c * v - 1.386294361119890618834464);
            NumericalScalar s(minp + r - w);
            // Quick acceptance steps
            // 2.609437912434100374600759 = 1 + log(5)
            if (s + 2.609437912434100374600759 >= 5.0 * z) return (p1 == minp) ? w / (maxp + w) : maxp / (maxp + w);
            NumericalScalar t(log(z));
            if (s > t) return (p1 == minp) ? w / (maxp + w) : maxp / (maxp + w);
            // Acceptance step
            if (r + sum * log(sum / (maxp + w)) >= t) return (p1 == minp) ? w / (maxp + w) : maxp / (maxp + w);
            }
        } // End Cheng
      if (maxp < 1.0) // Atkinson and Whittaker 1
        {
          NumericalScalar t(1.0 / (1.0 + sqrt(maxp * (1.0 - maxp) / (minp * (1.0 - minp)))));
          NumericalScalar tc(1.0 - t);
          NumericalScalar p(maxp * t / (maxp * t + minp * tc));
          for (;;)
            {
            NumericalScalar u(RandomGenerator::Generate());
            NumericalScalar e(-log(RandomGenerator::Generate()));
            if (u <= p)
              {
                NumericalScalar x(t * pow(u / p, 1.0 / minp));
                // Acceptation test
                if (e >= (1.0 - maxp) * log((1.0 - x) / tc)) return ((p1 == minp) ? x : 1.0 - x);
              }
            else
              {
                NumericalScalar x(1.0 - tc * pow((1.0 - u) / (1.0 - p), 1.0 / maxp));
                // Acceptation test
                if (e >= (1.0 - minp) * log(x / t)) return ((p1 == minp) ? x : 1.0 - x);
              }
            }
        } // End Atkinson and Whittaker 1
      // Remaining case, Atkinson and Whittaker 2
      NumericalScalar t((minp > 1) ? (1.0 - minp) / (maxp + 1.0 - minp) : 0.5);
      NumericalScalar tc(1.0 - t);
      NumericalScalar p(maxp * t / (maxp * t + minp * pow(tc, maxp)));
      for (;;)
        {
          NumericalScalar u(RandomGenerator::Generate());
          NumericalScalar e(-log(RandomGenerator::Generate()));
          if (u <= p)
            {
            NumericalScalar x(t * pow(u / p, 1.0 / minp));
            if (e >= (1.0 - maxp) * log(1.0 - x)) return ((p1 == minp) ? x : 1.0 - x);
            }
          else
            {
            NumericalScalar x(1.0 - tc * pow((1.0 - u) / (1.0 - p), 1.0 / maxp));
            if (e >= (1.0 - minp) * log(x / t)) return ((p1 == minp) ? x : 1.0 - x);
            }
        } // End Atkinson and Whittaker 2
      } // End of rBeta

      /*******************************************************************************************************/
      /* Normalized Gamma distribution, i.e. with a PDF equals to x ^ (k - 1) . exp(-x) / gamma(k) . (x > 0) */
      /*******************************************************************************************************/
      /* CDF */
      NumericalScalar DistFunc::pGamma(const NumericalScalar k, const NumericalScalar x, const Bool tail)
      {
#ifdef USE_DCDFLIB
      NumericalScalar ans;
      NumericalScalar qans;
      int ind(0);
      DCDFLIB::gamma_inc((double*)&k, (double*)&x, (double*)&ans, (double*)&qans, &ind);
      return ans;
#else
      return SpecFunc::GammaInc(k, x) / SpecFunc::Gamma(k);
#endif
      }
      /* CDF Inverse */
      NumericalScalar DistFunc::qGamma(const NumericalScalar k, const NumericalScalar x, const Bool tail)
      {
#ifdef USE_DCDFLIB
      NumericalScalar X;
      NumericalScalar X0(0.0);
      NumericalScalar p(x);
      NumericalScalar q(1.0 - x);
      int ierr;
      DCDFLIB::gamma_inc_inv((double*)&k, (double*)&X, (double*)&X0, (double*)&p, (double*)&q, &ierr);
      return X;
#else
      return SpecFunc::GammaIncInv(k, x * SpecFunc::Gamma(k));
#endif
      }
      /* Random number generation
       We use the algorithm described in:
       George Marsaglia and Wai Wan Tsang, "A Simple Method for Generating Gamma
         Variables": ACM Transactions on Mathematical Software, Vol. 26, No. 3,
         September 2000, Pages 363-372
       with a small optimization on the beta that appears in the squeezing function (1 + beta * x^4)*exp(-x^2/2).
       We also add the special treatment of the case k < 1
      */
      NumericalScalar DistFunc::rGamma(const NumericalScalar k)
      {
      // Special case k < 1.0
      NumericalScalar correction(1.0);
      NumericalScalar alpha(k);
      if (alpha < 1.0)
        {
          correction = pow(RandomGenerator::Generate(), 1.0 / alpha);
          alpha++;
        }
      NumericalScalar d(alpha - 0.3333333333333333333333333);
      NumericalScalar c(1.0 / sqrt(9.0 * d));
      NumericalScalar x;
      NumericalScalar v;
      for (;;)
        {
          do
            {
            x = rNormal();
            v = 1.0 + c * x;
            } while (v <= 0.0);
          v = v * v * v;
          NumericalScalar u(RandomGenerator::Generate());
          NumericalScalar x2(x * x);
          // Quick acceptation test
          // 0.03431688782875261396035499 is the numerical solution of the squeezing
          // problem described in 
          if (u < 1.0 - 0.03431688782875261396035499 * x2 * x2) return correction * d * v;
          // Acceptation test
          if (log(u) < 0.5 * x2 + d * (1.0 - v + log(v))) return correction * d * v;
        }
      } // End of rGamma

      /****************************/
      /* Kolmogorov distribution. */
      /****************************/
      /* CDF
       Adapted from:
       George Marsaglia and Wai Wan Tsang and Jingbo Wang (2003),
       "Evaluating Kolmogorov's distribution".
       Journal of Statistical Software, Volume 8, 2003, Issue 18.
       URL: http://www.jstatsoft.org/v08/i18/.
      */
      NumericalScalar DistFunc::pKolmogorov(const UnsignedLong n, const NumericalScalar x, const Bool tail)
      {
      // If x <= 0, we are at the left of the support
      if (x <= 0.0) return (tail ? 1.0 : 0.0);
      // If x >= 1, we are at the right of the support
      if (x >= 1.0) return (tail ? 0.0 : 1.0);
      const NumericalScalar s(x * x * n);
      // Less accurate (7 digits only) but much faster in the right tail
      if (s > 7.14 || (s > 3.76 && n > 99))
        {
          // Here, result is the complementary CDF
          const NumericalScalar result(2 * exp(-(2.000071 + 0.331 / sqrt(n) + 1.409 / n) * s));
          return (tail ? result : 1.0 - result);
        }
      // Now, the exact computation
      const UnsignedLong k(n * x + 1);
      const UnsignedLong m(2 * k - 1);
      SquareMatrix H(m);
      // Regular part of H
      const NumericalScalar factor(exp(SpecFunc::LnGamma(n + 1.0) / n - log(n)));
      for (UnsignedLong j = 0; j < m; ++j)
        {
          UnsignedLong k(1);
          NumericalScalar hij(factor);
          for (UnsignedLong i = j; i < m; ++i)
            {
            hij /= k;
            H(i, j) = hij;
            ++k;
            }
        }
      // Borders and upper diagonal of H
      const NumericalScalar h(k - n * x);
      NumericalScalar hi(h);
      for (UnsignedLong i = 1; i < m; ++i)
        {
          // Upper diagonal
          H(i - 1, i) = factor;
          // Left border
          H(i - 1, 0) *= (1.0 - hi);
          // Bottom border
          H(m - 1, m - i) *= (1.0 - hi);
          hi *= h;
        }
      // Left bottom corner, now hi = h^m
      H(m - 1, 0) *= (h > 0.5 ? 1.0 - 2.0 * hi + pow(2.0 * h - 1.0, m) : 1.0 - 2.0 * hi);
      // Here, result is the CDF
      const NumericalScalar result(H.power(n).operator()(k - 1, k - 1));
      return (tail ? 1.0 - result : result);
      }

      NumericalScalar DistFunc::pKolmogorovAsymptotic(const NumericalScalar x, const Bool tail)
      {
      // If x <= 0, we are at the left of the support
      // In fact, in double precision the numerical range starts slightly above 0.04
      if (x <= 0.04) return (tail ? 0.0 : 1.0);
      // Small x case
      if (x < 1.0)
        {
          NumericalScalar sum(0.0);
          const NumericalScalar z(-M_PI * M_PI / (8.0 * x * x));
          const NumericalScalar w(log(x));
          // All the terms starting at kmax are negligible in double precision
          const UnsignedLong kmax(0.5 * (1.0 + sqrt(1.0 + log(Precision) / z)) + 1.0);
          for (UnsignedLong k = 1; k < kmax; k += 2)
            {
            sum += exp(k * k * z - w);
            }
          return (tail ? 1.0 - sum * sqrt(2.0 * M_PI) : sum * sqrt(2.0 * M_PI));
        }
      else
        {
          const NumericalScalar z(-2.0 * x * x);
          NumericalScalar s(1.0);
          UnsignedLong k(1);
          NumericalScalar oldValue(1.0);
          NumericalScalar newValue(0.0);
          while(fabs(oldValue - newValue) > Precision)
            {
            oldValue = newValue;
            newValue += 2.0 * s * exp(z * k * k);
            s *= -1.0;
            ++k;
            }
          return (tail ? newValue : 1.0 - newValue);
        }
      }

      /************************************************************************************************************/
      /* Normalized NonCentralStudent distribution, i.e. with a PDF equals to (eq. 31.15 p.516 of the reference): */
      /* exp(-delta^2 / 2) * (nu / (nu + x^2)) ^ ((nu + 1) / 2) / (sqrt(nu * Pi) * Gamma(nu / 2) * SUM            */
      /* where SUM = sum_0^inf Gamma((nu + j + 1) / 2) . (x . delta . sqrt(2 / (nu + x . x))) ^ j / Gamma(j + 1)  */
      /* The summation is done starting at the even index k that maximize:                                        */
      /* Gamma((nu + j + 1) / 2) . (x . delta . sqrt(2 / (nu + x . x))) ^ j / Gamma(j + 1)                        */
      /* with respect to j. It can be shown using Stirling's formula that 2*[z/4] is a good approximation of k.   */
      /* Reference:                                                                                               */
      /* Norman L. Johnson, Samuel Kotz, N. Balakrishnan, "Continuous univariate distributions volume 2", second  */
      /* edition, 1995, Wiley Inter-Science                                                                       */
      /************************************************************************************************************/
      /* PDF */
      NumericalScalar DistFunc::dNonCentralStudent(const NumericalScalar nu, const NumericalScalar delta, const NumericalScalar x)
      {
      // Early exit for delta == 0, central Student PDF
      if (fabs(delta) < Precision) return exp(SpecFunc::LnGamma(0.5 * nu + 0.5) - SpecFunc::LnGamma(0.5 * nu) - 0.5 * log(M_PI * nu) - (0.5 * nu + 0.5) * log(1.0 + x * x / nu));
      // Case delta <> 0
      NumericalScalar halfnu(0.5 * nu);
      NumericalScalar nup1_2(halfnu + 0.5);
      NumericalScalar logConstant(-0.5 * delta * delta - SpecFunc::LnGamma(halfnu) - 0.5 * log(M_PI * nu));
      // Early exit for x == 0
      if (fabs(x) < Precision) return exp(logConstant + SpecFunc::LnGamma(nup1_2));
      // For x <> 0
      NumericalScalar w(1.0 / (nu + x * x));
      NumericalScalar value(exp(logConstant + nup1_2 * log(nu * w)));
      // Special treatment for very low value to avoid NaNs due to 0.Inf
      if (value < Precision)
        {
          if (x < 0.0) value = nu / x * (pNonCentralStudent(nu + 2.0, delta, x * sqrt(1.0 + 2.0 / nu)) - pNonCentralStudent(nu, delta, x));
          else value = -nu / x * (pNonCentralStudent(nu + 2.0, -delta, -x * sqrt(1.0 + 2.0 / nu)) - pNonCentralStudent(nu, -delta, -x));
          return (value <= 0.0 ? 0.0 : value);
        }
      NumericalScalar y(x * delta * sqrt(2.0 * w));
      NumericalScalar z(y * y);
      // Start at even index that maximize the coefficient in the sum
      NumericalScalar k(2.0 * floor(0.25 * y * y));
      // Loop forward and backward starting from k
      // Initialization
      NumericalScalar klogz(k * log(z));
      NumericalScalar twokp1(2.0 * k + 1.0);
      NumericalScalar lng2kp1(SpecFunc::LnGamma(twokp1));
      NumericalScalar lng2kp2(lng2kp1 + log(twokp1));
      NumericalScalar pkf(exp(SpecFunc::LnGamma(halfnu + k + 0.5) + klogz - lng2kp1));
      NumericalScalar pkb(pkf);
      NumericalScalar qkf(exp(SpecFunc::LnGamma(halfnu + k + 1.0) + klogz - lng2kp2) * y);
      NumericalScalar qkb(qkf);
      NumericalScalar sum(pkf + qkf);
      NumericalScalar zo4(0.25 * z);
      NumericalScalar izo4(1.0 / zo4);
      NumericalScalar dsf;
      NumericalScalar dsb;
      NumericalScalar error(1.0);
      NumericalScalar kf(k);
      NumericalScalar kb(k);
#define FORWARD_ITERATION                                   \
      pkf *= (nu + 2.0 * kf + 1.0) / ((2.0 * kf + 1.0) * (kf + 1.0)) * zo4; \
      qkf *= (nu + 2.0 * kf + 2.0) / ((2.0 * kf + 3.0) * (kf + 1.0)) * zo4; \
      dsf = pkf + qkf;                                \
      sum += dsf;
#define BACKWARD_ITERATION                                  \
      pkb *= kb * (2.0 * kb - 1.0) / (nu + 2.0 * kb - 1.0) * izo4;      \
      qkb *= (2.0 * kb + 1.0) * kb / (nu + 2.0 * kb) * izo4;            \
      dsb = pkb + qkb;                                \
      sum += dsb;

      // Here, i is an UnsignedLong as it is only a loop counter
      UnsignedLong i(1);

      UnsignedLong imax = static_cast<UnsignedLong>(k);
      if (imax > MaximumIteration) imax = MaximumIteration;
      while((error > Precision * (fabs(sum) + Precision)) && (i <= imax))
        {
          FORWARD_ITERATION;
          BACKWARD_ITERATION;
          error = fabs(dsf + dsb);
          ++kf;
          --kb;
          ++i;
        }
      // Do we have to perform further forward iterations?
      while ((error > Precision * (fabs(sum) + Precision)) && (i <= MaximumIteration))
        {
          FORWARD_ITERATION;
          ++kf;
          error = fabs(dsf);
          ++i;
        }
#undef FORWARD_ITERATION
#undef BACKWARD_ITERATION
      if (error > Precision * (fabs(sum) + Precision)) Log::Warn(OSS() << "Warning: in DistFunc::dNonCentralStudent(nu, delta, x), no convergence after " << i << " iterations. Error is " << error * sum << " value is " << value + 0.5 * sum << " for nu=" << nu << ", delta=" << delta << " and x=" << x);
      value *= sum;
      // Clip to [0,+inf[ in order to get rid of small rounding error
      return (value <= 0.0 ? 0.0 : value);
      }
      /* CDF
       We use the algorithm described in:
       Denise Benton, K. Krishnamoorthy, "Computing discrete mixtures of continuous distributions: noncentral chisquare, noncentral t
       and the distribution of the square of the sample multiple correlation coefficient",
       Computational Statistics & Data Analysis, 43 (2003) pp 249-267
      */
      NumericalScalar DistFunc::pNonCentralStudent(const NumericalScalar nu, const NumericalScalar delta, const NumericalScalar x, const Bool tail)
      {
      // Special case when |delta| << 1 
      if (fabs(delta) < Precision) return pStudent(nu, x, tail);
      // Call tnd routine for the other cases
      NumericalScalar t(x);
      NumericalScalar del(delta);
      NumericalScalar df(nu);
      if (tail)
        {
          t = -t;
          del = -del;
        }
      return TND_F77(&t, &df, &del);
      }
      /* Random number generation
       We use a transformation method based on Gamma and Normal transformation:
       If N is Normal(delta, 1) distributed and G is Gamma(nu / 2) distributed, 
       sqrt(2 * nu) * N / sqrt(G) is distributed according to NonCentralStudent(nu, delta)
      */
      NumericalScalar DistFunc::rNonCentralStudent(const NumericalScalar nu, const NumericalScalar delta)
      {
      NumericalScalar n(rNormal() + delta);
      NumericalScalar g(rGamma(0.5 * nu));
      return sqrt(0.5 * nu / g) * n;
      }

      /**************************************************************************************/
      /* Normalized Normal distribution, i.e. with a PDF equals to exp(-x^2/2) / sqrt(2.Pi) */
      /**************************************************************************************/
      /* CDF */
      NumericalScalar DistFunc::pNormal(const NumericalScalar x, const Bool tail)
      {
#ifdef USE_DCDFLIB
      NumericalScalar result;
      NumericalScalar ccum;
      DCDFLIB::cumnor((double*)&x, (double*)&result, (double*)&ccum);
      if (tail) return ccum;
      return result;
#else
      if (tail) return 0.5 - 0.5 * erf(x * M_SQRT1_2);
      return 0.5 + 0.5 * erf(x * M_SQRT1_2);
#endif
      }

      /* CDF inverse */
      NumericalScalar DistFunc::qNormal(const NumericalScalar x, const Bool tail)
      {
      if (x == 0.0) return -0.5 * sqrt(4.0 * (log(SpecFunc::ISQRT2PI) - SpecFunc::LogMinNumericalScalar));
      if (x == 1.0) return  0.5 * sqrt(4.0 * (log(SpecFunc::ISQRT2PI) - SpecFunc::LogMinNumericalScalar));
#ifdef USE_DCDFLIB
      NumericalScalar p(x);
      NumericalScalar q(1.0 - x);
      if (tail) return -DCDFLIB::dinvnr((double*)&p, (double*)&q);
      return DCDFLIB::dinvnr((double*)&p, (double*)&q);
#else
      if (tail) return -M_SQRT2l * SpecFunc::ErfInv(2.0 * x - 1.0);
      return M_SQRT2l * SpecFunc::ErfInv(2.0 * x - 1.0);
#endif
      }

      /* Random number generation
       We use the improved ziggurat method, see:
         Doornik, J.A. (2005), "An Improved Ziggurat Method to Generate Normal
         Random Samples", mimeo, Nuffield College, University of Oxford,
         and www.doornik.com/research.
      */
      NumericalScalar DistFunc::rNormal()
      {
      for (;;)
        {
          NumericalScalar u(2.0 * RandomGenerator::Generate() - 1.0);
          UnsignedLong index(RandomGenerator::IntegerGenerate(DistFunc::NumberOfBandNormalZigurrat));
          /* Are we in a rectangular band of the ziggurat? */
          if (fabs(u) < DistFunc::NormalZigguratRatio[index])
            {
            return u * NormalZigguratAbscissa[index + 1];
            }
          /* No, we are either on a wedge or in the upper tail of the Normal distribution */
          /* Are we in the bottom band? Sample from the tail of the Normal distribution */
          if (index == DistFunc::NumberOfBandNormalZigurrat - 1)
            {
            NumericalScalar x;
            NumericalScalar y;
            /* Marsaglia method */
            do
              {
                    x = log(RandomGenerator::Generate()) / DistFunc::NormalZigguratTail;
                    y = log(RandomGenerator::Generate());
              } while (-(y + y) < x * x);
            return (u > 0.0) ? x - DistFunc::NormalZigguratTail : DistFunc::NormalZigguratTail - x;
            }
          /* Are we in the wedges? Basic rejection method */
          NumericalScalar xI(NormalZigguratAbscissa[index]);
          NumericalScalar xIp1(NormalZigguratAbscissa[index + 1]);
          NumericalScalar x(u * xIp1);
          NumericalScalar pdfX(  exp(-0.5 * x * x));
          NumericalScalar pdfI(  exp(-0.5 * xI * xI));
          NumericalScalar pdfIp1(exp(-0.5 * xIp1 * xIp1));
          if (RandomGenerator::Generate() * (pdfI - pdfIp1) < pdfX - pdfIp1)
            {
                return x;
            }
        }
      }

      /**********************************************************************************/
      /* Poisson distribution, i.e. with a PDF equals to exp(-lambda) . lambda ^ k / k! */
      /**********************************************************************************/
      /* Random number generation
       For the small values of lambda, we use the method of inversion by sequential search described in:
       Luc Devroye, "Non-Uniform RandomVariate Generation", Springer-Verlag, 1986, available online at:
       http://cg.scs.carleton.ca/~luc/nonuniformrandomvariates.zip
       and with the important errata at:
       http://cg.scs.carleton.ca/~luc/errors.pdf
       For the large values of lambda, we use the ratio of uniform method described in:
       E. Stadlober, "The ratio of uniforms approach for generating discrete random variates". Journal of Computational and Applied Mathematics, vol. 31, no. 1, 1990, pp. 181-189.
      */
      NumericalScalar DistFunc::rPoisson(const NumericalScalar lambda)
      {
      NumericalScalar mu(floor(lambda));
      // Small case. The bound 6 is quite arbitrary, but must be < 80 to avoid overflow.
      if (mu < 6)
        {
          NumericalScalar x(0.0);
          NumericalScalar sum(exp(-lambda));
          NumericalScalar prod(sum);
          NumericalScalar u(RandomGenerator::Generate());
          for (;;)
            {
            if (u <= sum) return x;
            x++;
            prod *= lambda / x;
            sum += prod;
            }
        } // Small case
      // Large case
      NumericalScalar hatCenter(lambda + 0.5);
      NumericalScalar mode(floor(lambda));
      NumericalScalar logLambda(log(lambda));
      NumericalScalar pdfMode(mode * logLambda - SpecFunc::LnGamma(mode + 1.0));
      // 2.943035529371538572764190 = 8 / e
      // 0.898916162058898740826254 = 3 - 2 sqr(3 / e)
      NumericalScalar hatWidth(sqrt(2.943035529371538572764190 * (lambda + 0.5)) + 0.898916162058898740826254);
      NumericalScalar safetyBound(hatCenter + 6.0 * hatWidth);
      for (;;)
        {
          NumericalScalar u(RandomGenerator::Generate());
          NumericalScalar x(hatCenter + hatWidth * (RandomGenerator::Generate() - 0.5) / u);
        if (x < 0 || x >= safetyBound) continue;
        NumericalScalar k(floor(x));
        NumericalScalar logPdf(k * logLambda - SpecFunc::LnGamma(k + 1.0) - pdfMode);
        // Quick acceptance
        if (logPdf >= u * (4.0 - u) - 3.0) return k;
        // Quick rejection
        if (u * (u - logPdf) > 1.0) continue;
        // Acceptance
        if (2.0 * log(u) <= logPdf) return k;
      }
      }

      /********************************************************************************************************************************/
      /* Normalized Student distribution, i.e. with a PDF equals to (1 + x^2 / nu)^(-(1 + nu) / 2) / (sqrt(nu) . Beta(1 / 2, nu / 2)) */
      /********************************************************************************************************************************/
      /* CDF */
      NumericalScalar DistFunc::pStudent(const NumericalScalar nu, const NumericalScalar x, const Bool tail)
      {
      if (nu == 1.0) return 0.5 + atan(x) * M_1_PI;
#ifdef USE_DCDFLIB
      NumericalScalar cum;
      NumericalScalar ccum;
      DCDFLIB::cumt((double*)&x, (double*)&nu, (double*)&cum, (double*)&ccum);
      if (tail) return ccum;
      return cum;
#else
      return SpecFunc::BetaInc(0.5 * nu, 0.5, x);
#endif
      }
      /* CDF inverse */
      NumericalScalar DistFunc::qStudent(const NumericalScalar nu, const NumericalScalar x, const Bool tail)
      {
      if (nu == 1.0) return tan((x - 0.5) * M_PI);
#ifdef USE_DCDFLIB
      int which(2);
      NumericalScalar t;
      NumericalScalar q(1.0 - x);
      int status;
      NumericalScalar bound;
      DCDFLIB::cdft(&which, (double*)&x, (double*)&q, (double*)&t, (double*)&nu, &status, (double*)&bound);
      return t;
#else
      throw NotYetImplementedException(HERE);
#endif
      }
      /* Random number generation
       We use a transformation method based on Gamma and Normal transformation:
       If N is Normal(0, 1) distributed and G is Gamma(nu / 2) distributed, 
       sqrt(2 * nu) * N / sqrt(G) is distributed according to Student(nu)
      */
      NumericalScalar DistFunc::rStudent(const NumericalScalar nu)
      {
      NumericalScalar n(rNormal());
      NumericalScalar g(rGamma(0.5 * nu));
      return sqrt(0.5 * nu / g) * n;
      }
    } /* namespace Distribution */
  } /* namespace Uncertainty */
} /* namespace OpenTURNS */

Generated by  Doxygen 1.6.0   Back to index