blitz Version 1.0.2
Loading...
Searching...
No Matches
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
25namespace ranlib {
26
27template<typename T = double, typename IRNG = defaultIRNG,
28 typename stateTag = defaultState>
29class F {
30public:
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
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
85protected:
88};
89
90}
91
92#endif // BZ_RANDOM_F
Definition: F.h:29
T random()
Definition: F.h:57
void setDF(T _dfn, T _dfd)
Definition: F.h:46
void seed(IRNG_int s, IRNG_int r)
Definition: F.h:72
F(T numeratorDF, T denominatorDF, unsigned int i)
Definition: F.h:39
Gamma< T, IRNG, stateTag > dgamma
Definition: F.h:86
T T_numtype
Definition: F.h:31
T mindenom
Definition: F.h:87
F(T numeratorDF, T denominatorDF)
Definition: F.h:33
T dfd
Definition: F.h:87
Gamma< T, IRNG, stateTag > ngamma
Definition: F.h:86
T dfn
Definition: F.h:87
Definition: gamma.h:46
Definition: beta.h:50
sharedState defaultState
Definition: default.h:55
unsigned int IRNG_int
Definition: default.h:57
MersenneTwister defaultIRNG
Definition: default.h:120