forked from OSchip/llvm-project
Avoid huge main() functions and huge arrays. Patch from STL@microsoft.com
llvm-svn: 273354
This commit is contained in:
parent
7a580d0a72
commit
8d0d842e7b
|
@ -29,447 +29,491 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
D d(5, .75);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(30, .03125);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, .25);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.03);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 0);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 1);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == -inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
|
||||
void
|
||||
test6()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(400, 0.5);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs(skew - x_skew) < 0.01);
|
||||
assert(std::abs(kurtosis - x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test7()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 0.5);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs(skew - x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test8()
|
||||
{
|
||||
const int N = 100000;
|
||||
std::mt19937 gen1;
|
||||
std::mt19937 gen2;
|
||||
|
||||
std::binomial_distribution<> dist1(5, 0.1);
|
||||
std::binomial_distribution<unsigned> dist2(5, 0.1);
|
||||
|
||||
for(int i = 0; i < N; ++i)
|
||||
assert(dist1(gen1) == dist2(gen2));
|
||||
}
|
||||
|
||||
void
|
||||
test9()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 0.005);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
|
||||
void
|
||||
test10()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 0);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
|
||||
void
|
||||
test11()
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 1);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == -inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
D d(5, .75);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(30, .03125);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, .25);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.03);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 0);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 1);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == -inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(400, 0.5);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs(skew - x_skew) < 0.01);
|
||||
assert(std::abs(kurtosis - x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 0.5);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs(skew - x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
const int N = 100000;
|
||||
std::mt19937 gen1;
|
||||
std::mt19937 gen2;
|
||||
|
||||
std::binomial_distribution<> dist1(5, 0.1);
|
||||
std::binomial_distribution<unsigned> dist2(5, 0.1);
|
||||
|
||||
for(int i = 0; i < N; ++i)
|
||||
assert(dist1(gen1) == dist2(gen2));
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 0.005);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 0);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
{
|
||||
typedef std::binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0, 1);
|
||||
const int N = 100000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
// In this case:
|
||||
// skew computes to 0./0. == nan
|
||||
// kurtosis computes to 0./0. == nan
|
||||
// x_skew == -inf
|
||||
// x_kurtosis == inf
|
||||
// These tests are commented out because UBSan warns about division by 0
|
||||
// skew /= u.size() * dev * var;
|
||||
// kurtosis /= u.size() * var * var;
|
||||
// kurtosis -= 3;
|
||||
double x_mean = d.t() * d.p();
|
||||
double x_var = x_mean*(1-d.p());
|
||||
// double x_skew = (1-2*d.p()) / std::sqrt(x_var);
|
||||
// double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
// assert(skew == x_skew);
|
||||
// assert(kurtosis == x_kurtosis);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
test6();
|
||||
test7();
|
||||
test8();
|
||||
test9();
|
||||
test10();
|
||||
test11();
|
||||
}
|
||||
|
|
|
@ -29,246 +29,270 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(.03125);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.05);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::minstd_rand G;
|
||||
G g;
|
||||
D d(.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.75);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
|
||||
void
|
||||
test6()
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.96875);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(.03125);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.05);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::minstd_rand G;
|
||||
G g;
|
||||
D d(.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.75);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
{
|
||||
typedef std::geometric_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.96875);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
|
||||
double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
test6();
|
||||
}
|
||||
|
|
|
@ -29,244 +29,268 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::minstd_rand G;
|
||||
G g;
|
||||
D d(5, .25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(30, .03125);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, .25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 1);
|
||||
const int N = 1000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(400, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.04);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
|
||||
void
|
||||
test6()
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 0.05);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::minstd_rand G;
|
||||
G g;
|
||||
D d(5, .25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
|
||||
}
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(30, .03125);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, .25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(40, 1);
|
||||
const int N = 1000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(mean == x_mean);
|
||||
assert(var == x_var);
|
||||
}
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(400, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.04);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
{
|
||||
typedef std::negative_binomial_distribution<> D;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 0.05);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v <= d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(),
|
||||
double(0)) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.k() * (1 - d.p()) / d.p();
|
||||
double x_var = x_mean / d.p();
|
||||
double x_skew = (2 - d.p()) / std::sqrt(d.k() * (1 - d.p()));
|
||||
double x_kurtosis = 6. / d.k() + sqr(d.p()) / (d.k() * (1 - d.p()));
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
test6();
|
||||
}
|
||||
|
|
|
@ -29,216 +29,236 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./8192, 0.015625);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.05);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./32, 0.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./8, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.02);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.02);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.08);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.78125, 1.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.04);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.2);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./8192, 0.015625);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.05);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./32, 0.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-1./8, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.02);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.02);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.08);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.78125, 1.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(d.m() + sqr(d.s())/2);
|
||||
double x_var = (std::exp(sqr(d.s())) - 1) * std::exp(2*d.m() + sqr(d.s()));
|
||||
double x_skew = (std::exp(sqr(d.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(d.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(d.s())) + 2*std::exp(3*sqr(d.s())) +
|
||||
3*std::exp(2*sqr(d.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.04);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.2);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
}
|
||||
|
|
|
@ -29,222 +29,241 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./8192, 0.015625);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.05);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./32, 0.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./8, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.02);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(3, 4);
|
||||
P p;
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.02);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.08);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-0.78125, 1.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.04);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.2);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./8192, 0.015625);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.05);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.25);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./32, 0.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-1./8, 0.5);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.02);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.05);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(3, 4);
|
||||
P p;
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.02);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.08);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.4);
|
||||
}
|
||||
{
|
||||
typedef std::lognormal_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d;
|
||||
P p(-0.78125, 1.25);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
assert(v > 0);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = std::exp(p.m() + sqr(p.s())/2);
|
||||
double x_var = (std::exp(sqr(p.s())) - 1) * std::exp(2*p.m() + sqr(p.s()));
|
||||
double x_skew = (std::exp(sqr(p.s())) + 2) *
|
||||
std::sqrt((std::exp(sqr(p.s())) - 1));
|
||||
double x_kurtosis = std::exp(4*sqr(p.s())) + 2*std::exp(3*sqr(p.s())) +
|
||||
3*std::exp(2*sqr(p.s())) - 6;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.04);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.2);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.7);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
}
|
||||
|
|
|
@ -29,162 +29,178 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.5, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1.5, 3);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(3, 4);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(0.5, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(1.5, 3);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(3, 4);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = d.a() + d.b() * 0.577215665;
|
||||
double x_var = sqr(d.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
}
|
||||
|
|
|
@ -29,166 +29,182 @@ sqr(T x)
|
|||
return x * x;
|
||||
}
|
||||
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(0.5, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(1, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(1.5, 3);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(3, 4);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(0.5, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(1, 2);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(1.5, 3);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
{
|
||||
typedef std::extreme_value_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937 G;
|
||||
G g;
|
||||
D d(-0.5, 1);
|
||||
P p(3, 4);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g, p);
|
||||
u.push_back(v);
|
||||
}
|
||||
double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
|
||||
double var = 0;
|
||||
double skew = 0;
|
||||
double kurtosis = 0;
|
||||
for (int i = 0; i < u.size(); ++i)
|
||||
{
|
||||
double dbl = (u[i] - mean);
|
||||
double d2 = sqr(dbl);
|
||||
var += d2;
|
||||
skew += dbl * d2;
|
||||
kurtosis += d2 * d2;
|
||||
}
|
||||
var /= u.size();
|
||||
double dev = std::sqrt(var);
|
||||
skew /= u.size() * dev * var;
|
||||
kurtosis /= u.size() * var * var;
|
||||
kurtosis -= 3;
|
||||
double x_mean = p.a() + p.b() * 0.577215665;
|
||||
double x_var = sqr(p.b()) * 1.644934067;
|
||||
double x_skew = 1.139547;
|
||||
double x_kurtosis = 12./5;
|
||||
assert(std::abs((mean - x_mean) / x_mean) < 0.01);
|
||||
assert(std::abs((var - x_var) / x_var) < 0.01);
|
||||
assert(std::abs((skew - x_skew) / x_skew) < 0.01);
|
||||
assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
|
||||
}
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
}
|
||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -40,306 +40,330 @@ f(double x, double a, double m, double b, double c)
|
|||
return a + m*(sqr(x) - sqr(b))/2 + c*(x-b);
|
||||
}
|
||||
|
||||
int main()
|
||||
void
|
||||
test1()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {0, 1, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {0, 1, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {0, 0, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {1, 0, 0, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16};
|
||||
double p[] = {0, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14};
|
||||
double p[] = {1, 1};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {25, 62.5, 12.5, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
test2()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {0, 0, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
test3()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {1, 0, 0, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
test4()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16};
|
||||
double p[] = {0, 1, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
test5()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14};
|
||||
double p[] = {1, 1};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
test6()
|
||||
{
|
||||
typedef std::piecewise_linear_distribution<> D;
|
||||
typedef D::param_type P;
|
||||
typedef std::mt19937_64 G;
|
||||
G g;
|
||||
double b[] = {10, 14, 16, 17};
|
||||
double p[] = {25, 62.5, 12.5, 0};
|
||||
const size_t Np = sizeof(p) / sizeof(p[0]) - 1;
|
||||
D d(b, b+Np+1, p);
|
||||
const int N = 1000000;
|
||||
std::vector<D::result_type> u;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
D::result_type v = d(g);
|
||||
assert(d.min() <= v && v < d.max());
|
||||
u.push_back(v);
|
||||
}
|
||||
std::sort(u.begin(), u.end());
|
||||
int kp = -1;
|
||||
double a = std::numeric_limits<double>::quiet_NaN();
|
||||
double m = std::numeric_limits<double>::quiet_NaN();
|
||||
double bk = std::numeric_limits<double>::quiet_NaN();
|
||||
double c = std::numeric_limits<double>::quiet_NaN();
|
||||
std::vector<double> areas(Np);
|
||||
double S = 0;
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
{
|
||||
areas[i] = (p[i]+p[i+1])*(b[i+1]-b[i])/2;
|
||||
S += areas[i];
|
||||
}
|
||||
for (int i = 0; i < areas.size(); ++i)
|
||||
areas[i] /= S;
|
||||
for (int i = 0; i < Np+1; ++i)
|
||||
p[i] /= S;
|
||||
for (int i = 0; i < N; ++i)
|
||||
{
|
||||
int k = std::lower_bound(b, b+Np+1, u[i]) - b - 1;
|
||||
if (k != kp)
|
||||
{
|
||||
a = 0;
|
||||
for (int j = 0; j < k; ++j)
|
||||
a += areas[j];
|
||||
m = (p[k+1] - p[k]) / (b[k+1] - b[k]);
|
||||
bk = b[k];
|
||||
c = (b[k+1]*p[k] - b[k]*p[k+1]) / (b[k+1] - b[k]);
|
||||
kp = k;
|
||||
}
|
||||
assert(std::abs(f(u[i], a, m, bk, c) - double(i)/N) < .001);
|
||||
}
|
||||
}
|
||||
|
||||
int main()
|
||||
{
|
||||
test1();
|
||||
test2();
|
||||
test3();
|
||||
test4();
|
||||
test5();
|
||||
test6();
|
||||
}
|
||||
|
|
|
@ -15,7 +15,7 @@
|
|||
// the following definitions in this order:
|
||||
// 1) is_swappable, is_nothrow_swappable
|
||||
// 2) iter_swap, swap_ranges
|
||||
// 3) swap(T (&)[N], T(&)[N]
|
||||
// 3) swap(T (&)[N], T (&)[N])
|
||||
// This test checks that (1) and (2) see forward declarations
|
||||
// for (3).
|
||||
#include <type_traits>
|
||||
|
@ -27,7 +27,7 @@
|
|||
int main()
|
||||
{
|
||||
// Use a builtin type so we don't get ADL lookup.
|
||||
typedef double T[42][50];
|
||||
typedef double T[17][29];
|
||||
{
|
||||
LIBCPP_STATIC_ASSERT(std::__is_swappable<T>::value, "");
|
||||
#if TEST_STD_VER > 14
|
||||
|
@ -38,6 +38,6 @@ int main()
|
|||
T t1 = {};
|
||||
T t2 = {};
|
||||
std::iter_swap(t1, t2);
|
||||
std::swap_ranges(t1, t1 + 42, t2);
|
||||
std::swap_ranges(t1, t1 + 17, t2);
|
||||
}
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue