1
2 /* Web Polygraph http://www.web-polygraph.org/
3 * Copyright 2003-2011 The Measurement Factory
4 * Licensed under the Apache License, Version 2.0 */
5
6 #include "xstd/xstd.h"
7
8 #include "xstd/h/iostream.h"
9
10 #include "xstd/Assert.h"
11 #include "xstd/rndDistrs.h"
12
13
UnifDistr(RndGen * aGen,double aLo,double aHi)14 UnifDistr::UnifDistr(RndGen *aGen, double aLo, double aHi):
15 RndDistr(aGen), theLo(aLo), theHi(aHi) {
16
17 Assert(theHi >= theLo);
18 }
19
print(ostream & os,ArgPrinter p) const20 ostream &UnifDistr::print(ostream &os, ArgPrinter p) const {
21 os << pdfName() << '(';
22 p(os, theLo, 0);
23 p(os << ", ", theHi, 1);
24 return os << ')';
25 }
26
trial()27 double ExpDistr::trial() {
28 return -theMean * log(theGen->trial());
29 }
30
trial()31 double NormDistr::trial() {
32 double sum = 0;
33
34 for (int i = 0; i < 12; ++i)
35 sum += theGen->trial();
36
37 return theMean + theSDev*(sum - 6.0);
38 }
39
print(ostream & os,ArgPrinter p) const40 ostream &NormDistr::print(ostream &os, ArgPrinter p) const {
41 os << pdfName() << '(';
42 p(os, theMean, 0);
43 p(os << ", ", theSDev, 1);
44 return os << ')';
45 }
46
47 // translates mean and std dev into "internal" mu and sigma
48 // then creates the LognDistr object
ViaMean(RndGen * aGen,double aMean,double aSDev)49 LognDistr *LognDistr::ViaMean(RndGen *aGen, double aMean, double aSDev) {
50
51 const double twoLnM = 2 * log(aMean);
52 const double lnSum = log(aSDev*aSDev + aMean*aMean);
53
54 const double aMu = twoLnM - lnSum/2;
55 const double aSigmaSq = lnSum - twoLnM;
56
57 return new LognDistr(aGen, aMu, aSigmaSq, aMean, aSDev);
58 }
59
ViaMu(RndGen * aGen,double aMu,double aSigmaSq)60 LognDistr *LognDistr::ViaMu(RndGen *aGen, double aMu, double aSigmaSq) {
61 const double realMean = ::exp(aMu + aSigmaSq/2);
62 const double w = ::exp(aSigmaSq);
63 const double realSDev = sqrt(w * (w - 1) * ::exp(2*aMu)); // ?
64 return new LognDistr(aGen, aMu, aSigmaSq, realMean, realSDev);
65 }
66
LognDistr(RndGen * generator,double aMu,double aSigmaSq,double aMean,double aSDev)67 LognDistr::LognDistr(RndGen *generator, double aMu, double aSigmaSq, double aMean, double aSDev):
68 NormDistr(generator, aMu, sqrt(aSigmaSq)),
69 theRealMean(aMean), theRealSDev(aSDev) {
70 }
71
trial()72 double LognDistr::trial() {
73 return ::exp(NormDistr::trial());
74 }
75
76 // translates "internal" mu and sigma into mean
mean() const77 double LognDistr::mean() const {
78 return theRealMean;
79 }
80
81 // translates "internal" mu and sigma into std dev
sdev() const82 double LognDistr::sdev() const {
83 return theRealSDev;
84 }
85
print(ostream & os,ArgPrinter p) const86 ostream &LognDistr::print(ostream &os, ArgPrinter p) const {
87 os << pdfName() << '(';
88 p(os, theRealMean, 0);
89 p(os << ", ", theRealSDev, 1);
90 return os << ')';
91 }
92
93
94 /* Zipf */
95
96 // Euler-Mascheroni constant
97 static const double TheEMConst = 0.57721566490153286060651209;
98
trial()99 double ZipfDistr::trial() {
100 // return floor(pow(theWorldCap+1, theGen->trial()));
101 const double rn = theGen->trial();
102 return floor(pow(theWorldCap+1,rn));
103 }
104
ltrial(int min,int max)105 int ZipfDistr::ltrial(int min, int max) {
106 // wrong: (min, max]
107 // return min + (int) floor(pow(max-min+1, theGen->trial()));
108 // wrong: [min, max)
109 // [min, max]
110 min--;
111 const double rn = theGen->trial();
112 return min + (int) floor(pow(max-min+1, rn));
113 }
114
omega() const115 double ZipfDistr::omega() const {
116 // harmonic sum approximation: log(i) + gamma + 1/(2i)
117 return log((double)theWorldCap) + TheEMConst + 1./(2*theWorldCap);
118 }
119
print(ostream & os,ArgPrinter p) const120 ostream &ZipfDistr::print(ostream &os, ArgPrinter p) const {
121 os << pdfName() << '(';
122 p(os, theWorldCap, 0);
123 return os << ')';
124 }
125
126
127 /* Seq */
128
SeqDistr(RndGen * aGen,int aFirstId)129 SeqDistr::SeqDistr(RndGen *aGen, int aFirstId): RndDistr(aGen),
130 theFirstId(aFirstId), theLastId(aFirstId-1) {
131 }
132
trial()133 double SeqDistr::trial() {
134 return ++theLastId;
135 }
136
print(ostream & os,ArgPrinter p) const137 ostream &SeqDistr::print(ostream &os, ArgPrinter p) const {
138 os << pdfName() << '(';
139 p(os, theFirstId, 0);
140 return os << ')';
141 }
142