blitz Version 0.10
random/F.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // $Id: F.h,v 1.4 2011/03/25 22:41:17 julianc Exp $
00003 
00004 /*
00005  * F distribution
00006  *
00007  * This code has been adapted from RANDLIB.C 1.3, by
00008  * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
00009  * Code was originally by Ahrens and Dieter (see above).
00010  *
00011  * Adapter's notes:
00012  * BZ_NEEDS_WORK: how to handle seeding for the two gamma RNGs if
00013  * independentState is used?
00014  * BZ_NEEDS_WORK: code for handling possible overflow when xden
00015  * is tiny seems a bit flaky.
00016  */
00017 
00018 #ifndef BZ_RANDOM_F
00019 #define BZ_RANDOM_F
00020 
00021 #ifndef BZ_RANDOM_GAMMA
00022  #include <random/gamma.h>
00023 #endif
00024 
00025 BZ_NAMESPACE(ranlib)
00026 
00027 template<typename T = double, typename IRNG = defaultIRNG, 
00028     typename stateTag = defaultState>
00029 class F {
00030 public:
00031     typedef T T_numtype;
00032 
00033     F(T numeratorDF, T denominatorDF)
00034     {
00035         setDF(numeratorDF, denominatorDF);
00036         mindenom = 0.085 * blitz::tiny(T());
00037     }
00038 
00039   F(T numeratorDF, T denominatorDF, unsigned int i) :
00040     ngamma(i), dgamma(i)
00041     {
00042         setDF(numeratorDF, denominatorDF);
00043         mindenom = 0.085 * blitz::tiny(T());
00044     }
00045 
00046     void setDF(T _dfn, T _dfd)
00047     {
00048         BZPRECONDITION(_dfn > 0.0);
00049         BZPRECONDITION(_dfd > 0.0);
00050         dfn = _dfn;
00051         dfd = _dfd;
00052 
00053         ngamma.setMean(dfn/2.0);
00054         dgamma.setMean(dfd/2.0);
00055     }
00056 
00057     T random()
00058     {
00059         T xnum = 2.0 * ngamma.random() / dfn;
00060         T xden = 2.0 * ngamma.random() / dfd;
00061 
00062         // Rare event: Will an overflow probably occur?
00063         if (xden <= mindenom)
00064         {
00065             // Yes, just return huge(T())
00066             return blitz::huge(T());
00067         }
00068 
00069         return xnum / xden;
00070     }
00071 
00072   void seed(IRNG_int s, IRNG_int r)
00073     {
00074         // This is such a bad idea if independentState is used. Ugh.
00075         // If sharedState is used, it is merely inefficient (the
00076         // same RNG is seeded twice).
00077 
00078       // yes it's unacceptable -- changed to using two seeds / Patrik
00079       // in fact should probably be two uncorrelated IRNGs...
00080 
00081         ngamma.seed(s);
00082         dgamma.seed(r);
00083     }
00084 
00085 protected:
00086     Gamma<T,IRNG,stateTag> ngamma, dgamma;
00087     T dfn, dfd, mindenom;
00088 };
00089 
00090 BZ_NAMESPACE_END
00091 
00092 #endif // BZ_RANDOM_F
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines