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