blitz Version 0.10
|
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