00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include <it/random.h>
00027 #include <it/source.h>
00028 #include <it/source_func.h>
00029 #include <it/io.h>
00030 #include <it/math.h>
00031
00032
00033
00034 bvec source_binary (idx_t size, double p0)
00035 {
00036 bvec v = bvec_new (size);
00037 idx_t i;
00038
00039 for (i = 0; i < size; i++)
00040 v[i] = (it_rand () >= p0);
00041
00042 return (v);
00043 }
00044
00045
00046
00047 vec source_uniform (idx_t size, double a, double b)
00048 {
00049 vec v = vec_new (size);
00050 idx_t i;
00051
00052 for (i = 0; i < size; i++)
00053 v[i] = a + (b - a) * it_rand ();
00054
00055 return (v);
00056 }
00057
00058
00059 ivec source_uniform_int (idx_t size, int a, int b)
00060 {
00061 ivec v = ivec_new (size);
00062 idx_t i;
00063
00064 for (i = 0; i < size; i++)
00065 v[i] = floor ((double)a + (b - a) * it_rand ());
00066
00067 return (v);
00068 }
00069
00070
00071
00072 vec source_gaussian (idx_t size, double mean, double std)
00073 {
00074 vec v = vec_new (size);
00075 idx_t i;
00076
00077 for (i = 0; i < size; i++)
00078 v[i] = mean + std * it_randn ();
00079
00080 return (v);
00081 }
00082
00083
00084
00085
00086
00087
00088 vec source_pdf (idx_t size, double a, double b, it_function_t pdf,
00089 it_args_t args)
00090 {
00091 idx_t i;
00092 vec v = vec_new (size);
00093
00094 for (i = 0; i < size; i++)
00095 v[i] = it_randpdf (a, b, pdf, args);
00096
00097 return (v);
00098 }
00099
00100
00101
00102 ivec source_memoryless (idx_t size, vec pdf)
00103 {
00104 idx_t t;
00105 double x;
00106 int s, e, h;
00107 vec pdfcs;
00108 ivec S = ivec_new (size);
00109
00110 it_assert (is_valid_pdf (pdf, IT_EPSILON), "Invalid probability law");
00111 pdfcs = vec_cum_sum (pdf);
00112 vec_ins (pdfcs, 0, 0.0);
00113
00114 for (t = 0; t < size; t++) {
00115
00116
00117
00118
00119 x = it_rand ();
00120 s = 0;
00121 e = vec_length (pdfcs) - 1;
00122
00123 while (e > s + 1) {
00124
00125
00126 h = (s + e) / 2;
00127
00128 if (x < pdfcs[h])
00129 e = h;
00130 else
00131 s = h;
00132 }
00133
00134 if (x < pdfcs[e])
00135 S[t] = s;
00136 else
00137 S[t] = e;
00138 }
00139 vec_delete (pdfcs);
00140 return S;
00141 }