1 //===----------------------------------------------------------------------===//
2 //
3 // The LLVM Compiler Infrastructure
4 //
5 // This file is dual licensed under the MIT and the University of Illinois Open
6 // Source Licenses. See LICENSE.TXT for details.
7 //
8 //===----------------------------------------------------------------------===//
9 //
10 // REQUIRES: long_tests
11
12 // <random>
13
14 // template<class RealType = double>
15 // class piecewise_constant_distribution
16
17 // template<class _URNG> result_type operator()(_URNG& g);
18
19 #include <random>
20 #include <vector>
21 #include <iterator>
22 #include <numeric>
23 #include <algorithm> // for sort
24 #include <cassert>
25
26 template <class T>
27 inline
28 T
sqr(T x)29 sqr(T x)
30 {
31 return x*x;
32 }
33
34 void
test1()35 test1()
36 {
37 typedef std::piecewise_constant_distribution<> D;
38 typedef std::mt19937_64 G;
39 G g;
40 double b[] = {10, 14, 16, 17};
41 double p[] = {25, 62.5, 12.5};
42 const size_t Np = sizeof(p) / sizeof(p[0]);
43 D d(b, b+Np+1, p);
44 const int N = 1000000;
45 std::vector<D::result_type> u;
46 for (int i = 0; i < N; ++i)
47 {
48 D::result_type v = d(g);
49 assert(d.min() <= v && v < d.max());
50 u.push_back(v);
51 }
52 std::vector<double> prob(std::begin(p), std::end(p));
53 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
54 for (std::size_t i = 0; i < prob.size(); ++i)
55 prob[i] /= s;
56 std::sort(u.begin(), u.end());
57 for (std::size_t i = 0; i < Np; ++i)
58 {
59 typedef std::vector<D::result_type>::iterator I;
60 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
61 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
62 const size_t Ni = ub - lb;
63 if (prob[i] == 0)
64 assert(Ni == 0);
65 else
66 {
67 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
68 double mean = std::accumulate(lb, ub, 0.0) / Ni;
69 double var = 0;
70 double skew = 0;
71 double kurtosis = 0;
72 for (I j = lb; j != ub; ++j)
73 {
74 double dbl = (*j - mean);
75 double d2 = sqr(dbl);
76 var += d2;
77 skew += dbl * d2;
78 kurtosis += d2 * d2;
79 }
80 var /= Ni;
81 double dev = std::sqrt(var);
82 skew /= Ni * dev * var;
83 kurtosis /= Ni * var * var;
84 kurtosis -= 3;
85 double x_mean = (b[i+1] + b[i]) / 2;
86 double x_var = sqr(b[i+1] - b[i]) / 12;
87 double x_skew = 0;
88 double x_kurtosis = -6./5;
89 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
90 assert(std::abs((var - x_var) / x_var) < 0.01);
91 assert(std::abs(skew - x_skew) < 0.01);
92 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
93 }
94 }
95 }
96
97 void
test2()98 test2()
99 {
100 typedef std::piecewise_constant_distribution<> D;
101 typedef std::mt19937_64 G;
102 G g;
103 double b[] = {10, 14, 16, 17};
104 double p[] = {0, 62.5, 12.5};
105 const size_t Np = sizeof(p) / sizeof(p[0]);
106 D d(b, b+Np+1, p);
107 const int N = 1000000;
108 std::vector<D::result_type> u;
109 for (int i = 0; i < N; ++i)
110 {
111 D::result_type v = d(g);
112 assert(d.min() <= v && v < d.max());
113 u.push_back(v);
114 }
115 std::vector<double> prob(std::begin(p), std::end(p));
116 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
117 for (std::size_t i = 0; i < prob.size(); ++i)
118 prob[i] /= s;
119 std::sort(u.begin(), u.end());
120 for (std::size_t i = 0; i < Np; ++i)
121 {
122 typedef std::vector<D::result_type>::iterator I;
123 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
124 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
125 const size_t Ni = ub - lb;
126 if (prob[i] == 0)
127 assert(Ni == 0);
128 else
129 {
130 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
131 double mean = std::accumulate(lb, ub, 0.0) / Ni;
132 double var = 0;
133 double skew = 0;
134 double kurtosis = 0;
135 for (I j = lb; j != ub; ++j)
136 {
137 double dbl = (*j - mean);
138 double d2 = sqr(dbl);
139 var += d2;
140 skew += dbl * d2;
141 kurtosis += d2 * d2;
142 }
143 var /= Ni;
144 double dev = std::sqrt(var);
145 skew /= Ni * dev * var;
146 kurtosis /= Ni * var * var;
147 kurtosis -= 3;
148 double x_mean = (b[i+1] + b[i]) / 2;
149 double x_var = sqr(b[i+1] - b[i]) / 12;
150 double x_skew = 0;
151 double x_kurtosis = -6./5;
152 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
153 assert(std::abs((var - x_var) / x_var) < 0.01);
154 assert(std::abs(skew - x_skew) < 0.01);
155 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
156 }
157 }
158 }
159
160 void
test3()161 test3()
162 {
163 typedef std::piecewise_constant_distribution<> D;
164 typedef std::mt19937_64 G;
165 G g;
166 double b[] = {10, 14, 16, 17};
167 double p[] = {25, 0, 12.5};
168 const size_t Np = sizeof(p) / sizeof(p[0]);
169 D d(b, b+Np+1, p);
170 const int N = 1000000;
171 std::vector<D::result_type> u;
172 for (int i = 0; i < N; ++i)
173 {
174 D::result_type v = d(g);
175 assert(d.min() <= v && v < d.max());
176 u.push_back(v);
177 }
178 std::vector<double> prob(std::begin(p), std::end(p));
179 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
180 for (std::size_t i = 0; i < prob.size(); ++i)
181 prob[i] /= s;
182 std::sort(u.begin(), u.end());
183 for (std::size_t i = 0; i < Np; ++i)
184 {
185 typedef std::vector<D::result_type>::iterator I;
186 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
187 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
188 const size_t Ni = ub - lb;
189 if (prob[i] == 0)
190 assert(Ni == 0);
191 else
192 {
193 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
194 double mean = std::accumulate(lb, ub, 0.0) / Ni;
195 double var = 0;
196 double skew = 0;
197 double kurtosis = 0;
198 for (I j = lb; j != ub; ++j)
199 {
200 double dbl = (*j - mean);
201 double d2 = sqr(dbl);
202 var += d2;
203 skew += dbl * d2;
204 kurtosis += d2 * d2;
205 }
206 var /= Ni;
207 double dev = std::sqrt(var);
208 skew /= Ni * dev * var;
209 kurtosis /= Ni * var * var;
210 kurtosis -= 3;
211 double x_mean = (b[i+1] + b[i]) / 2;
212 double x_var = sqr(b[i+1] - b[i]) / 12;
213 double x_skew = 0;
214 double x_kurtosis = -6./5;
215 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
216 assert(std::abs((var - x_var) / x_var) < 0.01);
217 assert(std::abs(skew - x_skew) < 0.01);
218 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
219 }
220 }
221 }
222
223 void
test4()224 test4()
225 {
226 typedef std::piecewise_constant_distribution<> D;
227 typedef std::mt19937_64 G;
228 G g;
229 double b[] = {10, 14, 16, 17};
230 double p[] = {25, 62.5, 0};
231 const size_t Np = sizeof(p) / sizeof(p[0]);
232 D d(b, b+Np+1, p);
233 const int N = 1000000;
234 std::vector<D::result_type> u;
235 for (int i = 0; i < N; ++i)
236 {
237 D::result_type v = d(g);
238 assert(d.min() <= v && v < d.max());
239 u.push_back(v);
240 }
241 std::vector<double> prob(std::begin(p), std::end(p));
242 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
243 for (std::size_t i = 0; i < prob.size(); ++i)
244 prob[i] /= s;
245 std::sort(u.begin(), u.end());
246 for (std::size_t i = 0; i < Np; ++i)
247 {
248 typedef std::vector<D::result_type>::iterator I;
249 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
250 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
251 const size_t Ni = ub - lb;
252 if (prob[i] == 0)
253 assert(Ni == 0);
254 else
255 {
256 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
257 double mean = std::accumulate(lb, ub, 0.0) / Ni;
258 double var = 0;
259 double skew = 0;
260 double kurtosis = 0;
261 for (I j = lb; j != ub; ++j)
262 {
263 double dbl = (*j - mean);
264 double d2 = sqr(dbl);
265 var += d2;
266 skew += dbl * d2;
267 kurtosis += d2 * d2;
268 }
269 var /= Ni;
270 double dev = std::sqrt(var);
271 skew /= Ni * dev * var;
272 kurtosis /= Ni * var * var;
273 kurtosis -= 3;
274 double x_mean = (b[i+1] + b[i]) / 2;
275 double x_var = sqr(b[i+1] - b[i]) / 12;
276 double x_skew = 0;
277 double x_kurtosis = -6./5;
278 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
279 assert(std::abs((var - x_var) / x_var) < 0.01);
280 assert(std::abs(skew - x_skew) < 0.01);
281 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
282 }
283 }
284 }
285
286 void
test5()287 test5()
288 {
289 typedef std::piecewise_constant_distribution<> D;
290 typedef std::mt19937_64 G;
291 G g;
292 double b[] = {10, 14, 16, 17};
293 double p[] = {25, 0, 0};
294 const size_t Np = sizeof(p) / sizeof(p[0]);
295 D d(b, b+Np+1, p);
296 const int N = 100000;
297 std::vector<D::result_type> u;
298 for (int i = 0; i < N; ++i)
299 {
300 D::result_type v = d(g);
301 assert(d.min() <= v && v < d.max());
302 u.push_back(v);
303 }
304 std::vector<double> prob(std::begin(p), std::end(p));
305 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
306 for (std::size_t i = 0; i < prob.size(); ++i)
307 prob[i] /= s;
308 std::sort(u.begin(), u.end());
309 for (std::size_t i = 0; i < Np; ++i)
310 {
311 typedef std::vector<D::result_type>::iterator I;
312 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
313 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
314 const size_t Ni = ub - lb;
315 if (prob[i] == 0)
316 assert(Ni == 0);
317 else
318 {
319 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
320 double mean = std::accumulate(lb, ub, 0.0) / Ni;
321 double var = 0;
322 double skew = 0;
323 double kurtosis = 0;
324 for (I j = lb; j != ub; ++j)
325 {
326 double dbl = (*j - mean);
327 double d2 = sqr(dbl);
328 var += d2;
329 skew += dbl * d2;
330 kurtosis += d2 * d2;
331 }
332 var /= Ni;
333 double dev = std::sqrt(var);
334 skew /= Ni * dev * var;
335 kurtosis /= Ni * var * var;
336 kurtosis -= 3;
337 double x_mean = (b[i+1] + b[i]) / 2;
338 double x_var = sqr(b[i+1] - b[i]) / 12;
339 double x_skew = 0;
340 double x_kurtosis = -6./5;
341 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
342 assert(std::abs((var - x_var) / x_var) < 0.01);
343 assert(std::abs(skew - x_skew) < 0.01);
344 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
345 }
346 }
347 }
348
349 void
test6()350 test6()
351 {
352 typedef std::piecewise_constant_distribution<> D;
353 typedef std::mt19937_64 G;
354 G g;
355 double b[] = {10, 14, 16, 17};
356 double p[] = {0, 25, 0};
357 const size_t Np = sizeof(p) / sizeof(p[0]);
358 D d(b, b+Np+1, p);
359 const int N = 100000;
360 std::vector<D::result_type> u;
361 for (int i = 0; i < N; ++i)
362 {
363 D::result_type v = d(g);
364 assert(d.min() <= v && v < d.max());
365 u.push_back(v);
366 }
367 std::vector<double> prob(std::begin(p), std::end(p));
368 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
369 for (std::size_t i = 0; i < prob.size(); ++i)
370 prob[i] /= s;
371 std::sort(u.begin(), u.end());
372 for (std::size_t i = 0; i < Np; ++i)
373 {
374 typedef std::vector<D::result_type>::iterator I;
375 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
376 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
377 const size_t Ni = ub - lb;
378 if (prob[i] == 0)
379 assert(Ni == 0);
380 else
381 {
382 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
383 double mean = std::accumulate(lb, ub, 0.0) / Ni;
384 double var = 0;
385 double skew = 0;
386 double kurtosis = 0;
387 for (I j = lb; j != ub; ++j)
388 {
389 double dbl = (*j - mean);
390 double d2 = sqr(dbl);
391 var += d2;
392 skew += dbl * d2;
393 kurtosis += d2 * d2;
394 }
395 var /= Ni;
396 double dev = std::sqrt(var);
397 skew /= Ni * dev * var;
398 kurtosis /= Ni * var * var;
399 kurtosis -= 3;
400 double x_mean = (b[i+1] + b[i]) / 2;
401 double x_var = sqr(b[i+1] - b[i]) / 12;
402 double x_skew = 0;
403 double x_kurtosis = -6./5;
404 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
405 assert(std::abs((var - x_var) / x_var) < 0.01);
406 assert(std::abs(skew - x_skew) < 0.01);
407 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
408 }
409 }
410 }
411
412 void
test7()413 test7()
414 {
415 typedef std::piecewise_constant_distribution<> D;
416 typedef std::mt19937_64 G;
417 G g;
418 double b[] = {10, 14, 16, 17};
419 double p[] = {0, 0, 1};
420 const size_t Np = sizeof(p) / sizeof(p[0]);
421 D d(b, b+Np+1, p);
422 const int N = 100000;
423 std::vector<D::result_type> u;
424 for (int i = 0; i < N; ++i)
425 {
426 D::result_type v = d(g);
427 assert(d.min() <= v && v < d.max());
428 u.push_back(v);
429 }
430 std::vector<double> prob(std::begin(p), std::end(p));
431 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
432 for (std::size_t i = 0; i < prob.size(); ++i)
433 prob[i] /= s;
434 std::sort(u.begin(), u.end());
435 for (std::size_t i = 0; i < Np; ++i)
436 {
437 typedef std::vector<D::result_type>::iterator I;
438 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
439 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
440 const size_t Ni = ub - lb;
441 if (prob[i] == 0)
442 assert(Ni == 0);
443 else
444 {
445 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
446 double mean = std::accumulate(lb, ub, 0.0) / Ni;
447 double var = 0;
448 double skew = 0;
449 double kurtosis = 0;
450 for (I j = lb; j != ub; ++j)
451 {
452 double dbl = (*j - mean);
453 double d2 = sqr(dbl);
454 var += d2;
455 skew += dbl * d2;
456 kurtosis += d2 * d2;
457 }
458 var /= Ni;
459 double dev = std::sqrt(var);
460 skew /= Ni * dev * var;
461 kurtosis /= Ni * var * var;
462 kurtosis -= 3;
463 double x_mean = (b[i+1] + b[i]) / 2;
464 double x_var = sqr(b[i+1] - b[i]) / 12;
465 double x_skew = 0;
466 double x_kurtosis = -6./5;
467 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
468 assert(std::abs((var - x_var) / x_var) < 0.01);
469 assert(std::abs(skew - x_skew) < 0.01);
470 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
471 }
472 }
473 }
474
475 void
test8()476 test8()
477 {
478 typedef std::piecewise_constant_distribution<> D;
479 typedef std::mt19937_64 G;
480 G g;
481 double b[] = {10, 14, 16};
482 double p[] = {75, 25};
483 const size_t Np = sizeof(p) / sizeof(p[0]);
484 D d(b, b+Np+1, p);
485 const int N = 100000;
486 std::vector<D::result_type> u;
487 for (int i = 0; i < N; ++i)
488 {
489 D::result_type v = d(g);
490 assert(d.min() <= v && v < d.max());
491 u.push_back(v);
492 }
493 std::vector<double> prob(std::begin(p), std::end(p));
494 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
495 for (std::size_t i = 0; i < prob.size(); ++i)
496 prob[i] /= s;
497 std::sort(u.begin(), u.end());
498 for (std::size_t i = 0; i < Np; ++i)
499 {
500 typedef std::vector<D::result_type>::iterator I;
501 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
502 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
503 const size_t Ni = ub - lb;
504 if (prob[i] == 0)
505 assert(Ni == 0);
506 else
507 {
508 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
509 double mean = std::accumulate(lb, ub, 0.0) / Ni;
510 double var = 0;
511 double skew = 0;
512 double kurtosis = 0;
513 for (I j = lb; j != ub; ++j)
514 {
515 double dbl = (*j - mean);
516 double d2 = sqr(dbl);
517 var += d2;
518 skew += dbl * d2;
519 kurtosis += d2 * d2;
520 }
521 var /= Ni;
522 double dev = std::sqrt(var);
523 skew /= Ni * dev * var;
524 kurtosis /= Ni * var * var;
525 kurtosis -= 3;
526 double x_mean = (b[i+1] + b[i]) / 2;
527 double x_var = sqr(b[i+1] - b[i]) / 12;
528 double x_skew = 0;
529 double x_kurtosis = -6./5;
530 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
531 assert(std::abs((var - x_var) / x_var) < 0.01);
532 assert(std::abs(skew - x_skew) < 0.01);
533 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
534 }
535 }
536 }
537
538 void
test9()539 test9()
540 {
541 typedef std::piecewise_constant_distribution<> D;
542 typedef std::mt19937_64 G;
543 G g;
544 double b[] = {10, 14, 16};
545 double p[] = {0, 25};
546 const size_t Np = sizeof(p) / sizeof(p[0]);
547 D d(b, b+Np+1, p);
548 const int N = 100000;
549 std::vector<D::result_type> u;
550 for (int i = 0; i < N; ++i)
551 {
552 D::result_type v = d(g);
553 assert(d.min() <= v && v < d.max());
554 u.push_back(v);
555 }
556 std::vector<double> prob(std::begin(p), std::end(p));
557 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
558 for (std::size_t i = 0; i < prob.size(); ++i)
559 prob[i] /= s;
560 std::sort(u.begin(), u.end());
561 for (std::size_t i = 0; i < Np; ++i)
562 {
563 typedef std::vector<D::result_type>::iterator I;
564 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
565 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
566 const size_t Ni = ub - lb;
567 if (prob[i] == 0)
568 assert(Ni == 0);
569 else
570 {
571 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
572 double mean = std::accumulate(lb, ub, 0.0) / Ni;
573 double var = 0;
574 double skew = 0;
575 double kurtosis = 0;
576 for (I j = lb; j != ub; ++j)
577 {
578 double dbl = (*j - mean);
579 double d2 = sqr(dbl);
580 var += d2;
581 skew += dbl * d2;
582 kurtosis += d2 * d2;
583 }
584 var /= Ni;
585 double dev = std::sqrt(var);
586 skew /= Ni * dev * var;
587 kurtosis /= Ni * var * var;
588 kurtosis -= 3;
589 double x_mean = (b[i+1] + b[i]) / 2;
590 double x_var = sqr(b[i+1] - b[i]) / 12;
591 double x_skew = 0;
592 double x_kurtosis = -6./5;
593 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
594 assert(std::abs((var - x_var) / x_var) < 0.01);
595 assert(std::abs(skew - x_skew) < 0.01);
596 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
597 }
598 }
599 }
600
601 void
test10()602 test10()
603 {
604 typedef std::piecewise_constant_distribution<> D;
605 typedef std::mt19937_64 G;
606 G g;
607 double b[] = {10, 14, 16};
608 double p[] = {1, 0};
609 const size_t Np = sizeof(p) / sizeof(p[0]);
610 D d(b, b+Np+1, p);
611 const int N = 100000;
612 std::vector<D::result_type> u;
613 for (int i = 0; i < N; ++i)
614 {
615 D::result_type v = d(g);
616 assert(d.min() <= v && v < d.max());
617 u.push_back(v);
618 }
619 std::vector<double> prob(std::begin(p), std::end(p));
620 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
621 for (std::size_t i = 0; i < prob.size(); ++i)
622 prob[i] /= s;
623 std::sort(u.begin(), u.end());
624 for (std::size_t i = 0; i < Np; ++i)
625 {
626 typedef std::vector<D::result_type>::iterator I;
627 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
628 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
629 const size_t Ni = ub - lb;
630 if (prob[i] == 0)
631 assert(Ni == 0);
632 else
633 {
634 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
635 double mean = std::accumulate(lb, ub, 0.0) / Ni;
636 double var = 0;
637 double skew = 0;
638 double kurtosis = 0;
639 for (I j = lb; j != ub; ++j)
640 {
641 double dbl = (*j - mean);
642 double d2 = sqr(dbl);
643 var += d2;
644 skew += dbl * d2;
645 kurtosis += d2 * d2;
646 }
647 var /= Ni;
648 double dev = std::sqrt(var);
649 skew /= Ni * dev * var;
650 kurtosis /= Ni * var * var;
651 kurtosis -= 3;
652 double x_mean = (b[i+1] + b[i]) / 2;
653 double x_var = sqr(b[i+1] - b[i]) / 12;
654 double x_skew = 0;
655 double x_kurtosis = -6./5;
656 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
657 assert(std::abs((var - x_var) / x_var) < 0.01);
658 assert(std::abs(skew - x_skew) < 0.01);
659 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
660 }
661 }
662 }
663
664 void
test11()665 test11()
666 {
667 typedef std::piecewise_constant_distribution<> D;
668 typedef std::mt19937_64 G;
669 G g;
670 double b[] = {10, 14};
671 double p[] = {1};
672 const size_t Np = sizeof(p) / sizeof(p[0]);
673 D d(b, b+Np+1, p);
674 const int N = 100000;
675 std::vector<D::result_type> u;
676 for (int i = 0; i < N; ++i)
677 {
678 D::result_type v = d(g);
679 assert(d.min() <= v && v < d.max());
680 u.push_back(v);
681 }
682 std::vector<double> prob(std::begin(p), std::end(p));
683 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
684 for (std::size_t i = 0; i < prob.size(); ++i)
685 prob[i] /= s;
686 std::sort(u.begin(), u.end());
687 for (std::size_t i = 0; i < Np; ++i)
688 {
689 typedef std::vector<D::result_type>::iterator I;
690 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
691 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
692 const size_t Ni = ub - lb;
693 if (prob[i] == 0)
694 assert(Ni == 0);
695 else
696 {
697 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
698 double mean = std::accumulate(lb, ub, 0.0) / Ni;
699 double var = 0;
700 double skew = 0;
701 double kurtosis = 0;
702 for (I j = lb; j != ub; ++j)
703 {
704 double dbl = (*j - mean);
705 double d2 = sqr(dbl);
706 var += d2;
707 skew += dbl * d2;
708 kurtosis += d2 * d2;
709 }
710 var /= Ni;
711 double dev = std::sqrt(var);
712 skew /= Ni * dev * var;
713 kurtosis /= Ni * var * var;
714 kurtosis -= 3;
715 double x_mean = (b[i+1] + b[i]) / 2;
716 double x_var = sqr(b[i+1] - b[i]) / 12;
717 double x_skew = 0;
718 double x_kurtosis = -6./5;
719 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
720 assert(std::abs((var - x_var) / x_var) < 0.01);
721 assert(std::abs(skew - x_skew) < 0.01);
722 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
723 }
724 }
725 }
726
main()727 int main()
728 {
729 test1();
730 test2();
731 test3();
732 test4();
733 test5();
734 test6();
735 test7();
736 test8();
737 test9();
738 test10();
739 test11();
740 }
741