1 //
2 // random_histogram_example.c
3 //
4 // This example tests the random number generators for different
5 // distributions.
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <time.h>
12 #include <math.h>
13 #include <getopt.h>
14 #include "liquid.h"
15
16 #define OUTPUT_FILENAME "random_histogram_example.m"
17
18
19 // print usage/help message
usage()20 void usage()
21 {
22 printf("random_histogram_example [options]\n");
23 printf(" -h : print usage\n");
24 printf(" -N : number of trials\n");
25 printf(" -n : number of histogram bins\n");
26 printf(" -d : distribution: {uniform, normal, exp, weib, gamma, nak, rice}\n");
27 printf(" -u : u UNIFORM: lower edge\n");
28 printf(" -v : v UNIFORM: upper edge\n");
29 printf(" -e : eta NORMAL: mean\n");
30 printf(" -s : sigma NORMAL: standard deviation\n");
31 printf(" -l : lambda EXPONENTIAL: decay factor\n");
32 printf(" -a : alpha WEIBULL: shape\n");
33 printf(" -b : beta WEIBULL: spread\n");
34 printf(" -g : gamma WEIBULL: threshold\n");
35 printf(" -A : alpha GAMMA: shape\n");
36 printf(" -B : beta GAMMA: spread\n");
37 printf(" -m : m NAKAGAMI: shape\n");
38 printf(" -o : omega NAKAGAMI: spread\n");
39 printf(" -K : K RICE-K: spread\n");
40 printf(" -O : omega RICE-K: spread\n");
41 }
42
main(int argc,char * argv[])43 int main(int argc, char*argv[])
44 {
45 srand(time(NULL));
46 unsigned long int num_trials = 100000; // number of trials
47 unsigned int num_bins = 30;
48 enum {
49 UNIFORM=0,
50 NORMAL,
51 EXPONENTIAL,
52 WEIBULL,
53 GAMMA,
54 NAKAGAMIM,
55 RICEK
56 } distribution=NORMAL;
57
58 // distribution parameters
59 float u = 0.0f; // UNIFORM: lower edge
60 float v = 1.0f; // UNIFORM: upper edge
61 float eta = 0.0f; // NORMAL: mean
62 float sigma = 1.0f; // NORMAL: standard deviation
63 float lambda = 3.0f; // EXPONENTIAL: decay factor
64 float alphaw = 1.0f; // WEIBULL: shape
65 float betaw = 1.0f; // WEIBULL: spread
66 float gammaw = 1.0f; // WEIBULL: threshold
67 float alphag = 4.5f; // GAMMA: shape
68 float betag = 1.0f; // GAMMA: spread
69 float m = 4.5f; // NAKAGAMI: shape factor
70 float omeganak = 1.0f; // NAKAGMAI: spread factor
71 float K = 4.0f; // RICE-K: K-factor (shape)
72 float omegarice = 1.0f; // RICE-K: spread factor
73
74 int dopt;
75 while ((dopt = getopt(argc,argv,"hN:n:d:u:v:e:s:l:a:b:g:A:B:m:o:K:O:")) != EOF) {
76 switch (dopt) {
77 case 'h':
78 usage();
79 return 0;
80 case 'N': num_trials = atoi(optarg); break;
81 case 'n': num_bins = atoi(optarg); break;
82 case 'd':
83 if (strcmp(optarg,"uniform")==0) distribution = UNIFORM;
84 else if (strcmp(optarg,"normal")==0) distribution = NORMAL;
85 else if (strcmp(optarg,"exp")==0) distribution = EXPONENTIAL;
86 else if (strcmp(optarg,"weib")==0) distribution = WEIBULL;
87 else if (strcmp(optarg,"gamma")==0) distribution = GAMMA;
88 else if (strcmp(optarg,"nak")==0) distribution = NAKAGAMIM;
89 else if (strcmp(optarg,"rice")==0) distribution = RICEK;
90 else {
91 fprintf(stderr,"error: %s, unknown/unsupported distribution '%s'\n", argv[0], optarg);
92 exit(1);
93 }
94 case 'u': u = atof(optarg); break;
95 case 'v': v = atof(optarg); break;
96 case 'e': eta = atof(optarg); break;
97 case 's': sigma = atof(optarg); break;
98 case 'l': lambda = atof(optarg); break;
99 case 'a': alphaw = atof(optarg); break;
100 case 'b': betaw = atof(optarg); break;
101 case 'g': gammaw = atof(optarg); break;
102 case 'A': alphag = atof(optarg); break;
103 case 'B': betag = atof(optarg); break;
104 case 'm': m = atof(optarg); break;
105 case 'o': omeganak = atof(optarg); break;
106 case 'K': K = atof(optarg); break;
107 case 'O': omegarice = atof(optarg); break;
108 default:
109 exit(1);
110 }
111 }
112
113 // validate input
114 if (num_bins == 0) {
115 fprintf(stderr,"error: %s, number of bins must be greater than zero\n", argv[0]);
116 exit(1);
117 } else if (num_trials == 0) {
118 fprintf(stderr,"error: %s, number of trials must be greater than zero\n", argv[0]);
119 exit(1);
120 }
121
122 float xmin = 0.0f;
123 float xmax = 1.0f;
124
125 unsigned long int i;
126
127 // make a guess at the histogram range so we don't need to
128 // store all the generated random variables in a giant array.
129 if (distribution == UNIFORM) {
130 xmin = u - 0.08*(v-u); // lower edge less 8% range
131 xmax = v + 0.08*(v-u); // upper edge plus 8% range
132 } else if (distribution == NORMAL) {
133 xmin = eta - 4.0f*sigma;
134 xmax = eta + 4.0f*sigma;
135 } else if (distribution == EXPONENTIAL) {
136 xmin = 0.0f;
137 xmax = 7.0f / lambda;
138 } else if (distribution == WEIBULL) {
139 xmin = gammaw;
140 xmax = gammaw + betaw*powf( -logf(1e-3f), 1.0f/alphaw );
141 } else if (distribution == GAMMA) {
142 xmin = 0.0f;
143 xmax = 6.5 * betag + 2.0*alphag;
144 } else if (distribution == NAKAGAMIM) {
145 xmin = 0.0f;
146 xmax = 1.5f*( powf(omeganak, 0.8f) + 1.0f/m );
147 } else if (distribution == RICEK) {
148 xmin = 0.0f;
149 xmax = 3.0f*logf(omegarice+1.0f) + 1.5f/(K+1.0f);
150 } else {
151 fprintf(stderr, "error: %s, unknown/unsupported distribution\n", argv[0]);
152 exit(1);
153 }
154
155 //
156 //float xspan = xmax - xmin;
157 float bin_width = (xmax - xmin) / (num_bins);
158
159 // initialize histogram
160 unsigned int hist[num_bins];
161 for (i=0; i<num_bins; i++)
162 hist[i] = 0;
163
164 // generate random variables
165 float x = 0.0f;
166 float m1 = 0.0f; // first moment
167 float m2 = 0.0f; // second moment
168 for (i=0; i<num_trials; i++) {
169 switch (distribution) {
170 case UNIFORM: x = randuf(u,v); break;
171 case NORMAL: x = sigma*randnf() + eta; break;
172 case EXPONENTIAL: x = randexpf(lambda); break;
173 case WEIBULL: x = randweibf(alphaw,betaw,gammaw); break;
174 case GAMMA: x = randgammaf(alphag,betag); break;
175 case NAKAGAMIM: x = randnakmf(m,omeganak); break;
176 case RICEK: x = randricekf(K,omegarice); break;
177 default:
178 fprintf(stderr,"error: %s, unknown/unsupported distribution\n", argv[0]);
179 exit(1);
180 }
181
182 // compute bin index
183 unsigned int index;
184 float ihat = num_bins * (x - xmin) / (xmax - xmin);
185 if (ihat < 0.0f)
186 index = 0;
187 else
188 index = (unsigned int)ihat;
189
190 if (index >= num_bins)
191 index = num_bins-1;
192
193 hist[index]++;
194
195 // update statistics
196 m1 += x; // first moment
197 m2 += x*x; // second moment
198 }
199
200 //
201 m1 /= (float)num_trials;
202 m2 /= (float)num_trials;
203
204 // compute expected distribution
205 unsigned int num_steps = 100;
206 float xstep = (xmax - xmin) / (num_steps - 1);
207 float f[num_steps];
208 float F[num_steps];
209 for (i=0; i<num_steps; i++) {
210 x = xmin + i*xstep;
211 switch (distribution) {
212 case UNIFORM:
213 f[i] = randuf_pdf(x,u,v);
214 F[i] = randuf_cdf(x,u,v);
215 break;
216 case NORMAL:
217 f[i] = randnf_pdf(x,eta,sigma);
218 F[i] = randnf_cdf(x,eta,sigma);
219 break;
220 case EXPONENTIAL:
221 f[i] = randexpf_pdf(x,lambda);
222 F[i] = randexpf_cdf(x,lambda);
223 break;
224 case WEIBULL:
225 f[i] = randweibf_pdf(x,alphaw,betaw,gammaw);
226 F[i] = randweibf_cdf(x,alphaw,betaw,gammaw);
227 break;
228 case GAMMA:
229 f[i] = randgammaf_pdf(x,alphag,betag);
230 F[i] = randgammaf_cdf(x,alphag,betag);
231 break;
232 case NAKAGAMIM:
233 f[i] = randnakmf_pdf(x,m,omeganak);
234 F[i] = randnakmf_cdf(x,m,omeganak);
235 break;
236 case RICEK:
237 f[i] = randricekf_pdf(x,K,omegarice);
238 F[i] = randricekf_cdf(x,K,omegarice);
239 break;
240 default:
241 fprintf(stderr,"error: %s, unknown/unsupported distribution\n", argv[0]);
242 exit(1);
243 }
244 }
245
246 // print results to screen
247 // find max(hist)
248 unsigned int hist_max = 0;
249 for (i=0; i<num_bins; i++)
250 hist_max = hist[i] > hist_max ? hist[i] : hist_max;
251
252 printf("%8s : %6s [%6s]\n", "x", "count", "prob.");
253 for (i=0; i<num_bins; i++) {
254 printf("%8.2f : %6u [%6.4f]", xmin + i*bin_width, hist[i], (float)hist[i] / (float)num_trials);
255
256 unsigned int k;
257 unsigned int n = round(60 * (float)hist[i] / (float)hist_max);
258 for (k=0; k<n; k++)
259 printf("#");
260 printf("\n");
261 }
262
263 // print distribution info, statistics
264 printf("statistics:\n");
265 switch (distribution) {
266 case UNIFORM:
267 printf(" distribution : %s\n", "uniform");
268 printf(" u : %f\n", u);
269 printf(" v : %f\n", v);
270 break;
271 case NORMAL:
272 printf(" distribution : %s\n", "normal (Gauss)");
273 printf(" eta : %f\n", eta);
274 printf(" sigma : %f\n", sigma);
275 break;
276 case EXPONENTIAL:
277 printf(" distribution : %s\n", "exponential");
278 printf(" lambda : %f\n", lambda);
279 break;
280 case WEIBULL:
281 printf(" distribution : %s\n", "Weibull");
282 printf(" alpha : %f\n", alphaw);
283 printf(" beta : %f\n", betaw);
284 printf(" gamma : %f\n", gammaw);
285 break;
286 case GAMMA:
287 printf(" distribution : %s\n", "gamma");
288 printf(" alpha : %f\n", alphag);
289 printf(" beta : %f\n", betag);
290 break;
291 case NAKAGAMIM:
292 printf(" distribution : %s\n", "Nakagami-m");
293 printf(" m : %f\n", m);
294 printf(" omega : %f\n", omeganak);
295 break;
296 case RICEK:
297 printf(" distribution : %s\n", "Rice-K");
298 printf(" K : %f\n", K);
299 printf(" omega : %f\n", omegarice);
300 break;
301 default:
302 fprintf(stderr,"error: %s, unknown/unsupported distribution\n", argv[0]);
303 exit(1);
304 }
305 printf("\n");
306 printf(" samples : %8lu\n", num_trials);
307 printf(" first moment, E( x } : %8.3f\n", m1);
308 printf(" second moment, E{x^2} : %8.3f\n", m2);
309 printf(" variance : %8.3f\n", m2 - m1*m1);
310 printf(" standard deviation : %8.3f\n", sqrtf(m2 - m1*m1));
311
312 //
313 // export results
314 //
315 FILE * fid = fopen(OUTPUT_FILENAME,"w");
316 fprintf(fid,"%% %s : auto-generated file\n\n", OUTPUT_FILENAME);
317 fprintf(fid,"clear all;\n");
318 fprintf(fid,"close all;\n");
319 fprintf(fid,"xmin = %12.4e;\n", xmin);
320 fprintf(fid,"xmax = %12.4e;\n", xmax);
321 fprintf(fid,"num_bins = %u;\n", num_bins);
322 fprintf(fid,"xspan = xmax - xmin;\n");
323
324 float F_hat = 0.0f;
325 for (i=0; i<num_bins; i++) {
326 x = xmin + ((float)i + 0.5f)*bin_width;
327 float h = (float)(hist[i]) / (num_trials * bin_width);
328 fprintf(fid,"xh(%3lu) = %12.4e; h(%3lu) = %12.4e;\n", i+1, x, i+1, h);
329
330 x = xmin + ((float)i + 1.0f)*bin_width;
331 F_hat += h;
332 fprintf(fid,"xH(%3lu) = %12.4e; H(%3lu) = %12.4e;\n", i+1, x, i+1, F_hat);
333 }
334 fprintf(fid,"H = H/H(end);\n");
335
336 for (i=0; i<num_steps; i++) {
337 x = xmin + i*xstep;
338 fprintf(fid,"xf(%3lu) = %12.4e; f(%3lu) = %12.4e; F(%3lu) = %12.4e;\n", i+1, x, i+1, f[i], i+1, F[i]);
339 }
340
341 // plot results
342 fprintf(fid,"figure;\n");
343 fprintf(fid,"plot(xh,h,'x', xf,f,'-');\n");
344 fprintf(fid,"xlabel('x');\n");
345 fprintf(fid,"ylabel('f_x(x)');\n");
346 fprintf(fid,"axis([(xmin-0.1*xspan) (xmax+0.1*xspan) 0 1.1*max([h f])]);\n");
347 fprintf(fid,"legend('histogram','true PDF',1);\n");
348
349 // plot results
350 fprintf(fid,"figure;\n");
351 fprintf(fid,"plot(xH,H,'x', xf,F,'-');\n");
352 fprintf(fid,"xlabel('x');\n");
353 fprintf(fid,"ylabel('f_x(x)');\n");
354 //fprintf(fid,"axis([(xmin-0.1*xspan) (xmax+0.1*xspan) 0 1]);\n");
355 fprintf(fid,"legend('histogram','true CDF',0);\n");
356
357 fclose(fid);
358 printf("results written to %s.\n",OUTPUT_FILENAME);
359
360
361 printf("done.\n");
362 return 0;
363 }
364
365