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