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

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.samp/rand.dist.samp.pconst/eval.pass.cpp