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