blitz  Version 1.0.2
F.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 // $Id$
3 
4 /*
5  * F distribution
6  *
7  * This code has been adapted from RANDLIB.C 1.3, by
8  * Barry W. Brown, James Lovato, Kathy Russell, and John Venier.
9  * Code was originally by Ahrens and Dieter (see above).
10  *
11  * Adapter's notes:
12  * BZ_NEEDS_WORK: how to handle seeding for the two gamma RNGs if
13  * independentState is used?
14  * BZ_NEEDS_WORK: code for handling possible overflow when xden
15  * is tiny seems a bit flaky.
16  */
17 
18 #ifndef BZ_RANDOM_F
19 #define BZ_RANDOM_F
20 
21 #ifndef BZ_RANDOM_GAMMA
22  #include <random/gamma.h>
23 #endif
24 
25 namespace ranlib {
26 
27 template<typename T = double, typename IRNG = defaultIRNG,
28  typename stateTag = defaultState>
29 class F {
30 public:
31  typedef T T_numtype;
32 
33  F(T numeratorDF, T denominatorDF)
34  {
35  setDF(numeratorDF, denominatorDF);
36  mindenom = 0.085 * blitz::tiny(T());
37  }
38 
39  F(T numeratorDF, T denominatorDF, unsigned int i) :
40  ngamma(i), dgamma(i)
41  {
42  setDF(numeratorDF, denominatorDF);
43  mindenom = 0.085 * blitz::tiny(T());
44  }
45 
46  void setDF(T _dfn, T _dfd)
47  {
48  BZPRECONDITION(_dfn > 0.0);
49  BZPRECONDITION(_dfd > 0.0);
50  dfn = _dfn;
51  dfd = _dfd;
52 
53  ngamma.setMean(dfn/2.0);
54  dgamma.setMean(dfd/2.0);
55  }
56 
57  T random()
58  {
59  T xnum = 2.0 * ngamma.random() / dfn;
60  T xden = 2.0 * ngamma.random() / dfd;
61 
62  // Rare event: Will an overflow probably occur?
63  if (xden <= mindenom)
64  {
65  // Yes, just return huge(T())
66  return blitz::huge(T());
67  }
68 
69  return xnum / xden;
70  }
71 
73  {
74  // This is such a bad idea if independentState is used. Ugh.
75  // If sharedState is used, it is merely inefficient (the
76  // same RNG is seeded twice).
77 
78  // yes it's unacceptable -- changed to using two seeds / Patrik
79  // in fact should probably be two uncorrelated IRNGs...
80 
81  ngamma.seed(s);
82  dgamma.seed(r);
83  }
84 
85 protected:
88 };
89 
90 }
91 
92 #endif // BZ_RANDOM_F
_bz_global blitz::IndexPlaceholder< 9 > r
Definition: indexexpr.h:265
Definition: beta.h:50
Definition: gamma.h:45
_bz_global blitz::IndexPlaceholder< 0 > i
Definition: indexexpr.h:256
T dfd
Definition: F.h:87
Gamma< T, IRNG, stateTag > ngamma
Definition: F.h:86
T T_numtype
Definition: F.h:31
T dfn
Definition: F.h:87
Definition: F.h:29
F(T numeratorDF, T denominatorDF)
Definition: F.h:33
void setDF(T _dfn, T _dfd)
Definition: F.h:46
F(T numeratorDF, T denominatorDF, unsigned int i)
Definition: F.h:39
unsigned int IRNG_int
Definition: default.h:57
MersenneTwister defaultIRNG
Definition: default.h:120
T random()
Definition: F.h:57
T mindenom
Definition: F.h:87
void seed(IRNG_int s, IRNG_int r)
Definition: F.h:72
_bz_global blitz::IndexPlaceholder< 10 > s
Definition: indexexpr.h:266
sharedState defaultState
Definition: default.h:55
Gamma< T, IRNG, stateTag > dgamma
Definition: F.h:86