|
2928
|
1 /* |
|
|
2 |
|
|
3 Copyright (C) 1996, 1997 John W. Eaton |
|
|
4 |
|
|
5 This file is part of Octave. |
|
|
6 |
|
|
7 Octave is free software; you can redistribute it and/or modify it |
|
|
8 under the terms of the GNU General Public License as published by the |
|
|
9 Free Software Foundation; either version 2, or (at your option) any |
|
|
10 later version. |
|
|
11 |
|
|
12 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
15 for more details. |
|
|
16 |
|
|
17 You should have received a copy of the GNU General Public License |
|
|
18 along with Octave; see the file COPYING. If not, write to the Free |
|
5307
|
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
20 02110-1301, USA. |
|
2928
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include <config.h> |
|
|
26 #endif |
|
|
27 |
|
|
28 #include <ctime> |
|
|
29 |
|
|
30 #include <string> |
|
|
31 |
|
|
32 #include "f77-fcn.h" |
|
|
33 #include "lo-mappers.h" |
|
4307
|
34 #include "oct-rand.h" |
|
4153
|
35 #include "quit.h" |
|
2928
|
36 |
|
|
37 #include "defun-dld.h" |
|
|
38 #include "error.h" |
|
|
39 #include "gripes.h" |
|
|
40 #include "oct-obj.h" |
|
|
41 #include "unwind-prot.h" |
|
|
42 #include "utils.h" |
|
|
43 |
|
4307
|
44 static octave_value |
|
5730
|
45 do_rand (const octave_value_list& args, int nargin, const char *fcn, |
|
|
46 bool additional_arg = false) |
|
2928
|
47 { |
|
4307
|
48 octave_value retval; |
|
5730
|
49 NDArray a; |
|
|
50 int idx = 0; |
|
|
51 dim_vector dims; |
|
2928
|
52 |
|
5730
|
53 if (additional_arg) |
|
|
54 { |
|
|
55 if (nargin == 0) |
|
|
56 { |
|
|
57 error ("%s: expecting at least one argument", fcn); |
|
|
58 goto done; |
|
|
59 } |
|
|
60 else if (args(0).is_string()) |
|
|
61 additional_arg = false; |
|
|
62 else |
|
|
63 { |
|
|
64 a = args(0).array_value (); |
|
|
65 if (error_state) |
|
|
66 { |
|
|
67 error ("%s: expecting scalar or matrix arguments", fcn); |
|
|
68 goto done; |
|
|
69 } |
|
|
70 idx++; |
|
|
71 nargin--; |
|
|
72 } |
|
|
73 } |
|
2928
|
74 |
|
4543
|
75 switch (nargin) |
|
2928
|
76 { |
|
4543
|
77 case 0: |
|
|
78 { |
|
5730
|
79 if (additional_arg) |
|
|
80 dims = a.dims (); |
|
|
81 else |
|
|
82 { |
|
|
83 dims.resize (2); |
|
4543
|
84 |
|
5730
|
85 dims(0) = 1; |
|
|
86 dims(1) = 1; |
|
|
87 } |
|
4543
|
88 goto gen_matrix; |
|
|
89 } |
|
|
90 break; |
|
2928
|
91 |
|
4543
|
92 case 1: |
|
|
93 { |
|
5730
|
94 octave_value tmp = args(idx); |
|
4543
|
95 |
|
|
96 if (tmp.is_string ()) |
|
|
97 { |
|
|
98 std::string s_arg = tmp.string_value (); |
|
2928
|
99 |
|
4543
|
100 if (s_arg == "dist") |
|
|
101 { |
|
|
102 retval = octave_rand::distribution (); |
|
|
103 } |
|
5730
|
104 else if (s_arg == "seed") |
|
4543
|
105 { |
|
|
106 retval = octave_rand::seed (); |
|
|
107 } |
|
5730
|
108 else if (s_arg == "state" || s_arg == "twister") |
|
|
109 { |
|
|
110 retval = octave_rand::state (); |
|
|
111 } |
|
4543
|
112 else if (s_arg == "uniform") |
|
|
113 { |
|
|
114 octave_rand::uniform_distribution (); |
|
|
115 } |
|
|
116 else if (s_arg == "normal") |
|
|
117 { |
|
|
118 octave_rand::normal_distribution (); |
|
|
119 } |
|
5730
|
120 else if (s_arg == "exponential") |
|
|
121 { |
|
|
122 octave_rand::exponential_distribution (); |
|
|
123 } |
|
|
124 else if (s_arg == "poisson") |
|
|
125 { |
|
|
126 octave_rand::poisson_distribution (); |
|
|
127 } |
|
|
128 else if (s_arg == "gamma") |
|
|
129 { |
|
|
130 octave_rand::gamma_distribution (); |
|
|
131 } |
|
4543
|
132 else |
|
4664
|
133 error ("%s: unrecognized string argument", fcn); |
|
4543
|
134 } |
|
|
135 else if (tmp.is_scalar_type ()) |
|
|
136 { |
|
|
137 double dval = tmp.double_value (); |
|
2928
|
138 |
|
4543
|
139 if (xisnan (dval)) |
|
|
140 { |
|
4664
|
141 error ("%s: NaN is invalid a matrix dimension", fcn); |
|
4543
|
142 } |
|
|
143 else |
|
|
144 { |
|
|
145 dims.resize (2); |
|
|
146 |
|
5275
|
147 dims(0) = NINTbig (tmp.double_value ()); |
|
|
148 dims(1) = NINTbig (tmp.double_value ()); |
|
2928
|
149 |
|
4543
|
150 if (! error_state) |
|
|
151 goto gen_matrix; |
|
|
152 } |
|
|
153 } |
|
|
154 else if (tmp.is_range ()) |
|
|
155 { |
|
|
156 Range r = tmp.range_value (); |
|
|
157 |
|
|
158 if (r.all_elements_are_ints ()) |
|
|
159 { |
|
5275
|
160 octave_idx_type n = r.nelem (); |
|
4543
|
161 |
|
|
162 dims.resize (n); |
|
|
163 |
|
5275
|
164 octave_idx_type base = NINTbig (r.base ()); |
|
|
165 octave_idx_type incr = NINTbig (r.inc ()); |
|
|
166 octave_idx_type lim = NINTbig (r.limit ()); |
|
2928
|
167 |
|
4543
|
168 if (base < 0 || lim < 0) |
|
4664
|
169 error ("%s: all dimensions must be nonnegative", fcn); |
|
4543
|
170 else |
|
|
171 { |
|
5275
|
172 for (octave_idx_type i = 0; i < n; i++) |
|
4543
|
173 { |
|
|
174 dims(i) = base; |
|
|
175 base += incr; |
|
|
176 } |
|
2928
|
177 |
|
4543
|
178 goto gen_matrix; |
|
|
179 } |
|
|
180 } |
|
|
181 else |
|
4664
|
182 error ("%s: expecting all elements of range to be integers", |
|
|
183 fcn); |
|
4543
|
184 } |
|
|
185 else if (tmp.is_matrix_type ()) |
|
|
186 { |
|
|
187 Array<int> iv = tmp.int_vector_value (true); |
|
|
188 |
|
|
189 if (! error_state) |
|
|
190 { |
|
5275
|
191 octave_idx_type len = iv.length (); |
|
2928
|
192 |
|
4543
|
193 dims.resize (len); |
|
|
194 |
|
5275
|
195 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
196 { |
|
5275
|
197 octave_idx_type elt = iv(i); |
|
4543
|
198 |
|
|
199 if (elt < 0) |
|
|
200 { |
|
4664
|
201 error ("%s: all dimensions must be nonnegative", fcn); |
|
4543
|
202 goto done; |
|
|
203 } |
|
|
204 |
|
|
205 dims(i) = iv(i); |
|
|
206 } |
|
2928
|
207 |
|
4543
|
208 goto gen_matrix; |
|
|
209 } |
|
|
210 else |
|
4664
|
211 error ("%s: expecting integer vector", fcn); |
|
4543
|
212 } |
|
|
213 else |
|
|
214 { |
|
|
215 gripe_wrong_type_arg ("rand", tmp); |
|
|
216 return retval; |
|
|
217 } |
|
|
218 } |
|
|
219 break; |
|
|
220 |
|
|
221 default: |
|
|
222 { |
|
5730
|
223 octave_value tmp = args(idx); |
|
4543
|
224 |
|
|
225 if (nargin == 2 && tmp.is_string ()) |
|
|
226 { |
|
5164
|
227 std::string ts = tmp.string_value (); |
|
|
228 |
|
5730
|
229 if (ts == "seed") |
|
4543
|
230 { |
|
5782
|
231 if (args(idx+1).is_real_scalar ()) |
|
|
232 { |
|
|
233 double d = args(idx+1).double_value (); |
|
2928
|
234 |
|
5782
|
235 if (! error_state) |
|
|
236 octave_rand::seed (d); |
|
|
237 } |
|
|
238 else |
|
|
239 error ("%s: seed must be a real scalar", fcn); |
|
4543
|
240 } |
|
5730
|
241 else if (ts == "state" || ts == "twister") |
|
|
242 { |
|
|
243 ColumnVector s = |
|
|
244 ColumnVector (args(idx+1).vector_value(false, true)); |
|
|
245 |
|
|
246 if (! error_state) |
|
|
247 octave_rand::state (s); |
|
|
248 } |
|
4543
|
249 else |
|
4664
|
250 error ("%s: unrecognized string argument", fcn); |
|
4543
|
251 } |
|
|
252 else |
|
|
253 { |
|
|
254 dims.resize (nargin); |
|
|
255 |
|
|
256 for (int i = 0; i < nargin; i++) |
|
|
257 { |
|
5760
|
258 dims(i) = args(idx+i).int_value (); |
|
4543
|
259 |
|
|
260 if (error_state) |
|
|
261 { |
|
4664
|
262 error ("%s: expecting integer arguments", fcn); |
|
4543
|
263 goto done; |
|
|
264 } |
|
|
265 } |
|
|
266 |
|
|
267 goto gen_matrix; |
|
|
268 } |
|
|
269 } |
|
|
270 break; |
|
2928
|
271 } |
|
|
272 |
|
4543
|
273 done: |
|
2928
|
274 |
|
|
275 return retval; |
|
|
276 |
|
|
277 gen_matrix: |
|
|
278 |
|
5355
|
279 dims.chop_trailing_singletons (); |
|
|
280 |
|
5730
|
281 if (additional_arg) |
|
|
282 { |
|
|
283 if (a.length() == 1) |
|
|
284 return octave_rand::nd_array (dims, a(0)); |
|
|
285 else |
|
|
286 { |
|
|
287 if (a.dims() != dims) |
|
|
288 { |
|
|
289 error ("%s: mismatch in argument size", fcn); |
|
|
290 return retval; |
|
|
291 } |
|
|
292 octave_idx_type len = a.length (); |
|
|
293 NDArray m (dims); |
|
|
294 double *v = m.fortran_vec (); |
|
|
295 for (octave_idx_type i = 0; i < len; i++) |
|
|
296 v[i] = octave_rand::scalar (a(i)); |
|
|
297 return m; |
|
|
298 } |
|
|
299 } |
|
|
300 else |
|
|
301 return octave_rand::nd_array (dims); |
|
2928
|
302 } |
|
|
303 |
|
4665
|
304 DEFUN_DLD (rand, args, , |
|
3369
|
305 "-*- texinfo -*-\n\ |
|
|
306 @deftypefn {Loadable Function} {} rand (@var{x})\n\ |
|
|
307 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ |
|
5730
|
308 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ |
|
|
309 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ |
|
3369
|
310 Return a matrix with random elements uniformly distributed on the\n\ |
|
|
311 interval (0, 1). The arguments are handled the same as the arguments\n\ |
|
5730
|
312 for @code{eye}.\n\ |
|
|
313 \n\ |
|
|
314 You can query the state of the random number generator using the\n\ |
|
3369
|
315 form\n\ |
|
2928
|
316 \n\ |
|
3369
|
317 @example\n\ |
|
5730
|
318 v = rand (\"state\")\n\ |
|
|
319 @end example\n\ |
|
|
320 \n\ |
|
|
321 This returns a column vector @var{v} of length 625. Later, you can\n\ |
|
|
322 restore the random number generator to the state @var{v}\n\ |
|
|
323 using the form\n\ |
|
|
324 \n\ |
|
|
325 @example\n\ |
|
|
326 rand (\"state\", v)\n\ |
|
3369
|
327 @end example\n\ |
|
|
328 \n\ |
|
|
329 @noindent\n\ |
|
5730
|
330 You may also initialize the state vector from an arbitrary vector of\n\ |
|
|
331 length <= 625 for @var{v}. This new state will be a hash based on the\n\ |
|
|
332 the value of @var{v}, not @var{v} itself.\n\ |
|
|
333 \n\ |
|
|
334 By default, the generator is initialized from @code{/dev/urandom} if it is\n\ |
|
|
335 available, otherwise from cpu time, wall clock time and the current\n\ |
|
|
336 fraction of a second.\n\ |
|
|
337 \n\ |
|
|
338 @code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\ |
|
|
339 (See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ |
|
|
340 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ |
|
|
341 Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\ |
|
|
342 @url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\ |
|
|
343 Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\ |
|
|
344 values together, otherwise the generator state can be learned after\n\ |
|
|
345 reading 624 consecutive values.\n\ |
|
|
346 \n\ |
|
|
347 @code{rand} includes a second random number generator, that was the\n\ |
|
|
348 previous generator used in octave. The new generator is used by default\n\ |
|
|
349 as it is significantly faster than the old generator, and produces\n\ |
|
|
350 random numebrs with a significantly longer cycle time. However, in\n\ |
|
|
351 some circumstances it might be desireable to obtain the same random\n\ |
|
|
352 sequences as used by the old generators. To do this the keyword\n\ |
|
|
353 \"seed\" is used to specify that the old generators should be use,\n\ |
|
|
354 as in\n\ |
|
2928
|
355 \n\ |
|
3369
|
356 @example\n\ |
|
5730
|
357 rand (\"seed\", val)\n\ |
|
3369
|
358 @end example\n\ |
|
|
359 \n\ |
|
5730
|
360 which sets the seed of the generator to @var{val}. The seed of the\n\ |
|
|
361 generator can be queried with\n\ |
|
|
362 \n\ |
|
|
363 @example\n\ |
|
|
364 s = rand (\"seed\")\n\ |
|
|
365 @end example\n\ |
|
|
366 \n\ |
|
|
367 However, it should be noted that querying the seed will not cause\n\ |
|
|
368 @code{rand} to use the old generators, only setting the seed will.\n\ |
|
|
369 To cause @code{rand} to once again use the new generators, the\n\ |
|
|
370 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ |
|
|
371 @seealso{randn,rande,randg,randp}\n\ |
|
3369
|
372 @end deftypefn") |
|
2928
|
373 { |
|
4307
|
374 octave_value retval; |
|
2928
|
375 |
|
|
376 int nargin = args.length (); |
|
|
377 |
|
4543
|
378 retval = do_rand (args, nargin, "rand"); |
|
2928
|
379 |
|
|
380 return retval; |
|
|
381 } |
|
|
382 |
|
5730
|
383 /* |
|
|
384 %!test # 'state' can be a scalar |
|
|
385 %! rand('state',12); x = rand(1,4); |
|
|
386 %! rand('state',12); y = rand(1,4); |
|
|
387 %! assert(x,y); |
|
|
388 %!test # 'state' can be a vector |
|
|
389 %! rand('state',[12,13]); x=rand(1,4); |
|
|
390 %! rand('state',[12;13]); y=rand(1,4); |
|
|
391 %! assert(x,y); |
|
|
392 %!test # querying 'state' doesn't disturb sequence |
|
|
393 %! rand('state',12); rand(1,2); x=rand(1,2); |
|
|
394 %! rand('state',12); rand(1,2); |
|
|
395 %! s=rand('state'); y=rand(1,2); |
|
|
396 %! assert(x,y); |
|
|
397 %! rand('state',s); z=rand(1,2); |
|
|
398 %! assert(x,z); |
|
|
399 %!test # 'seed' must be a scalar |
|
|
400 %! rand('seed',12); x = rand(1,4); |
|
|
401 %! rand('seed',12); y = rand(1,4); |
|
|
402 %! assert(x,y); |
|
|
403 %!error(rand('seed',[12,13])) |
|
|
404 %!test # querying 'seed' returns a value which can be used later |
|
|
405 %! s=rand('seed'); x=rand(1,2); |
|
|
406 %! rand('seed',s); y=rand(1,2); |
|
|
407 %! assert(x,y); |
|
|
408 %!test # querying 'seed' doesn't disturb sequence |
|
|
409 %! rand('seed',12); rand(1,2); x=rand(1,2); |
|
|
410 %! rand('seed',12); rand(1,2); |
|
|
411 %! s=rand('seed'); y=rand(1,2); |
|
|
412 %! assert(x,y); |
|
|
413 %! rand('seed',s); z=rand(1,2); |
|
|
414 %! assert(x,z); |
|
|
415 */ |
|
|
416 |
|
|
417 /* |
|
|
418 %!test |
|
|
419 %! % statistical tests may fail occasionally. |
|
|
420 %! rand("state",12); |
|
|
421 %! x = rand(100000,1); |
|
|
422 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
423 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
424 %! assert(mean(x),0.5,0.0024); |
|
|
425 %! assert(var(x),1/48,0.0632); |
|
|
426 %! assert(skewness(x),0,0.012); |
|
|
427 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
428 %!test |
|
|
429 %! % statistical tests may fail occasionally. |
|
|
430 %! rand("seed",12); |
|
|
431 %! x = rand(100000,1); |
|
|
432 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
433 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
434 %! assert(mean(x),0.5,0.0024); |
|
|
435 %! assert(var(x),1/48,0.0632); |
|
|
436 %! assert(skewness(x),0,0.012); |
|
|
437 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
438 */ |
|
|
439 |
|
|
440 |
|
4307
|
441 static std::string current_distribution = octave_rand::distribution (); |
|
|
442 |
|
2928
|
443 static void |
|
|
444 reset_rand_generator (void *) |
|
|
445 { |
|
4307
|
446 octave_rand::distribution (current_distribution); |
|
2928
|
447 } |
|
|
448 |
|
4665
|
449 DEFUN_DLD (randn, args, , |
|
3369
|
450 "-*- texinfo -*-\n\ |
|
|
451 @deftypefn {Loadable Function} {} randn (@var{x})\n\ |
|
|
452 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ |
|
5730
|
453 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
|
|
454 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ |
|
|
455 Return a matrix with normally distributed random elements. The\n\ |
|
|
456 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
3369
|
457 \n\ |
|
5730
|
458 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
459 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ |
|
|
460 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
461 J. Statistical Software, vol 5, 2000,\n\ |
|
|
462 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
2928
|
463 \n\ |
|
5730
|
464 @seealso{rand,rande,randg,randp}\n\ |
|
3369
|
465 @end deftypefn") |
|
2928
|
466 { |
|
4307
|
467 octave_value retval; |
|
2928
|
468 |
|
|
469 int nargin = args.length (); |
|
|
470 |
|
4543
|
471 unwind_protect::begin_frame ("randn"); |
|
2928
|
472 |
|
4543
|
473 // This relies on the fact that elements are popped from the unwind |
|
|
474 // stack in the reverse of the order they are pushed |
|
|
475 // (i.e. current_distribution will be reset before calling |
|
|
476 // reset_rand_generator()). |
|
2928
|
477 |
|
4543
|
478 unwind_protect::add (reset_rand_generator, 0); |
|
|
479 unwind_protect_str (current_distribution); |
|
2928
|
480 |
|
4543
|
481 current_distribution = "normal"; |
|
2928
|
482 |
|
4543
|
483 octave_rand::distribution (current_distribution); |
|
2928
|
484 |
|
4543
|
485 retval = do_rand (args, nargin, "randn"); |
|
2928
|
486 |
|
4543
|
487 unwind_protect::run_frame ("randn"); |
|
2928
|
488 |
|
|
489 return retval; |
|
|
490 } |
|
|
491 |
|
|
492 /* |
|
5730
|
493 %!test |
|
|
494 %! % statistical tests may fail occasionally. |
|
|
495 %! rand("state",12); |
|
|
496 %! x = randn(100000,1); |
|
|
497 %! assert(mean(x),0,0.01); |
|
|
498 %! assert(var(x),1,0.02); |
|
|
499 %! assert(skewness(x),0,0.02); |
|
|
500 %! assert(kurtosis(x),0,0.04); |
|
|
501 %!test |
|
|
502 %! % statistical tests may fail occasionally. |
|
|
503 %! rand("seed",12); |
|
|
504 %! x = randn(100000,1); |
|
|
505 %! assert(mean(x),0,0.01); |
|
|
506 %! assert(var(x),1,0.02); |
|
|
507 %! assert(skewness(x),0,0.02); |
|
|
508 %! assert(kurtosis(x),0,0.04); |
|
|
509 */ |
|
|
510 |
|
|
511 DEFUN_DLD (rande, args, , |
|
|
512 "-*- texinfo -*-\n\ |
|
|
513 @deftypefn {Loadable Function} {} rande (@var{x})\n\ |
|
|
514 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ |
|
|
515 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ |
|
|
516 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ |
|
|
517 Return a matrix with exponentially distributed random elements. The\n\ |
|
|
518 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
|
519 \n\ |
|
|
520 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
521 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ |
|
|
522 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
523 J. Statistical Software, vol 5, 2000,\n\ |
|
|
524 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
|
525 @seealso{rand,randn,randg,randp}\n\ |
|
|
526 @end deftypefn") |
|
|
527 { |
|
|
528 octave_value retval; |
|
|
529 |
|
|
530 int nargin = args.length (); |
|
|
531 |
|
|
532 unwind_protect::begin_frame ("rande"); |
|
|
533 |
|
|
534 // This relies on the fact that elements are popped from the unwind |
|
|
535 // stack in the reverse of the order they are pushed |
|
|
536 // (i.e. current_distribution will be reset before calling |
|
|
537 // reset_rand_generator()). |
|
|
538 |
|
|
539 unwind_protect::add (reset_rand_generator, 0); |
|
|
540 unwind_protect_str (current_distribution); |
|
|
541 |
|
|
542 current_distribution = "exponential"; |
|
|
543 |
|
|
544 octave_rand::distribution (current_distribution); |
|
|
545 |
|
|
546 retval = do_rand (args, nargin, "rande"); |
|
|
547 |
|
|
548 unwind_protect::run_frame ("rande"); |
|
|
549 |
|
|
550 return retval; |
|
|
551 } |
|
|
552 |
|
|
553 /* |
|
|
554 %!test |
|
|
555 %! % statistical tests may fail occasionally |
|
|
556 %! rand("state",12); |
|
|
557 %! x = rande(100000,1); |
|
|
558 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
559 %! assert(mean(x),1,0.01); |
|
|
560 %! assert(var(x),1,0.03); |
|
|
561 %! assert(skewness(x),2,0.06); |
|
|
562 %! assert(kurtosis(x),6,0.7); |
|
|
563 %!test |
|
|
564 %! % statistical tests may fail occasionally |
|
|
565 %! rand("seed",12); |
|
|
566 %! x = rande(100000,1); |
|
|
567 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
568 %! assert(mean(x),1,0.01); |
|
|
569 %! assert(var(x),1,0.03); |
|
|
570 %! assert(skewness(x),2,0.06); |
|
|
571 %! assert(kurtosis(x),6,0.7); |
|
|
572 */ |
|
|
573 |
|
|
574 DEFUN_DLD (randg, args, , |
|
|
575 "-*- texinfo -*-\n\ |
|
|
576 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ |
|
|
577 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ |
|
|
578 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ |
|
|
579 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ |
|
|
580 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ |
|
|
581 The arguments are handled the same as the arguments for @code{rand},\n\ |
|
|
582 except for the argument @var{a}.\n\ |
|
|
583 \n\ |
|
|
584 This can be used to generate many distributions:\n\ |
|
|
585 \n\ |
|
|
586 @table @asis\n\ |
|
|
587 @item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\ |
|
|
588 @example\n\ |
|
|
589 r = b*randg(a)\n\ |
|
|
590 @end example\n\ |
|
|
591 @item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\ |
|
|
592 @example\n\ |
|
|
593 r1 = randg(a,1)\n\ |
|
|
594 r = r1 / (r1 + randg(b,1))\n\ |
|
|
595 @end example\n\ |
|
|
596 @item @code{Erlang(a, n)}\n\ |
|
|
597 @example\n\ |
|
|
598 r = a*randg(n)\n\ |
|
|
599 @end example\n\ |
|
|
600 @item @code{chisq(df)} for @code{df > 0}\n\ |
|
|
601 @example\n\ |
|
|
602 r = 2*randg(df/2)\n\ |
|
|
603 @end example\n\ |
|
|
604 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ |
|
|
605 @example\n\ |
|
|
606 r = randn() / sqrt(2*randg(df/2)/df)\n\ |
|
|
607 @end example\n\ |
|
|
608 @item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\ |
|
|
609 @example\n\ |
|
|
610 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\ |
|
|
611 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\ |
|
|
612 r = r1 / r2\n\n\ |
|
|
613 @end example\n\ |
|
|
614 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ |
|
|
615 @example\n\ |
|
|
616 r = randp((1-p)/p * randg(n))\n\ |
|
|
617 @end example\n\ |
|
|
618 @item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\ |
|
|
619 (use chisq if @code{L = 0})\n\ |
|
|
620 @example\n\ |
|
|
621 r = randp(L/2)\n\ |
|
|
622 r(r > 0) = 2*randg(r(r > 0))\n\ |
|
|
623 r(df > 0) += 2*randg(df(df > 0)/2)\n\ |
|
|
624 @end example\n\ |
|
|
625 @item @code{Dirichlet(a1,...,ak)}\n\ |
|
|
626 @example\n\ |
|
|
627 r = (randg(a1),...,randg(ak))\n\ |
|
|
628 r = r / sum(r)\n\ |
|
|
629 @end example\n\ |
|
|
630 @end table\n\ |
|
|
631 @seealso{rand,randn,rande,randp}\n\ |
|
|
632 @end deftypefn") |
|
|
633 { |
|
|
634 octave_value retval; |
|
|
635 |
|
|
636 int nargin = args.length (); |
|
|
637 |
|
|
638 if (nargin < 1) |
|
|
639 error ("randg: insufficient arguments"); |
|
|
640 else |
|
|
641 { |
|
|
642 unwind_protect::begin_frame ("randg"); |
|
|
643 |
|
|
644 // This relies on the fact that elements are popped from the unwind |
|
|
645 // stack in the reverse of the order they are pushed |
|
|
646 // (i.e. current_distribution will be reset before calling |
|
|
647 // reset_rand_generator()). |
|
|
648 |
|
|
649 unwind_protect::add (reset_rand_generator, 0); |
|
|
650 unwind_protect_str (current_distribution); |
|
|
651 |
|
|
652 current_distribution = "gamma"; |
|
|
653 |
|
|
654 octave_rand::distribution (current_distribution); |
|
|
655 |
|
|
656 retval = do_rand (args, nargin, "randg", true); |
|
|
657 |
|
|
658 unwind_protect::run_frame ("randg"); |
|
|
659 } |
|
|
660 |
|
|
661 return retval; |
|
|
662 } |
|
|
663 |
|
|
664 /* |
|
|
665 %!test |
|
|
666 %! rand("state",12) |
|
|
667 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
668 %!test |
|
|
669 %! % statistical tests may fail occasionally. |
|
|
670 %! a=0.1; x = randg(a,100000,1); |
|
|
671 %! assert(mean(x), a, 0.01); |
|
|
672 %! assert(var(x), a, 0.01); |
|
|
673 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
674 %! assert(kurtosis(x),6/a, 50.); |
|
|
675 %!test |
|
|
676 %! % statistical tests may fail occasionally. |
|
|
677 %! a=0.95; x = randg(a,100000,1); |
|
|
678 %! assert(mean(x), a, 0.01); |
|
|
679 %! assert(var(x), a, 0.04); |
|
|
680 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
681 %! assert(kurtosis(x),6/a, 2.); |
|
|
682 %!test |
|
|
683 %! % statistical tests may fail occasionally. |
|
|
684 %! a=1; x = randg(a,100000,1); |
|
|
685 %! assert(mean(x),a, 0.01); |
|
|
686 %! assert(var(x),a, 0.04); |
|
|
687 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
688 %! assert(kurtosis(x),6/a, 2.); |
|
|
689 %!test |
|
|
690 %! % statistical tests may fail occasionally. |
|
|
691 %! a=10; x = randg(a,100000,1); |
|
|
692 %! assert(mean(x), a, 0.1); |
|
|
693 %! assert(var(x), a, 0.5); |
|
|
694 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
695 %! assert(kurtosis(x),6/a, 0.5); |
|
|
696 %!test |
|
|
697 %! % statistical tests may fail occasionally. |
|
|
698 %! a=100; x = randg(a,100000,1); |
|
|
699 %! assert(mean(x), a, 0.2); |
|
|
700 %! assert(var(x), a, 2.); |
|
|
701 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
702 %! assert(kurtosis(x),6/a, 0.2); |
|
|
703 %!test |
|
|
704 %! rand("seed",12) |
|
|
705 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
706 %!test |
|
|
707 %! % statistical tests may fail occasionally. |
|
|
708 %! a=0.1; x = randg(a,100000,1); |
|
|
709 %! assert(mean(x), a, 0.01); |
|
|
710 %! assert(var(x), a, 0.01); |
|
|
711 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
712 %! assert(kurtosis(x),6/a, 50.); |
|
|
713 %!test |
|
|
714 %! % statistical tests may fail occasionally. |
|
|
715 %! a=0.95; x = randg(a,100000,1); |
|
|
716 %! assert(mean(x), a, 0.01); |
|
|
717 %! assert(var(x), a, 0.04); |
|
|
718 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
719 %! assert(kurtosis(x),6/a, 2.); |
|
|
720 %!test |
|
|
721 %! % statistical tests may fail occasionally. |
|
|
722 %! a=1; x = randg(a,100000,1); |
|
|
723 %! assert(mean(x),a, 0.01); |
|
|
724 %! assert(var(x),a, 0.04); |
|
|
725 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
726 %! assert(kurtosis(x),6/a, 2.); |
|
|
727 %!test |
|
|
728 %! % statistical tests may fail occasionally. |
|
|
729 %! a=10; x = randg(a,100000,1); |
|
|
730 %! assert(mean(x), a, 0.1); |
|
|
731 %! assert(var(x), a, 0.5); |
|
|
732 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
733 %! assert(kurtosis(x),6/a, 0.5); |
|
|
734 %!test |
|
|
735 %! % statistical tests may fail occasionally. |
|
|
736 %! a=100; x = randg(a,100000,1); |
|
|
737 %! assert(mean(x), a, 0.2); |
|
|
738 %! assert(var(x), a, 2.); |
|
|
739 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
740 %! assert(kurtosis(x),6/a, 0.2); |
|
|
741 */ |
|
|
742 |
|
|
743 |
|
|
744 DEFUN_DLD (randp, args, , |
|
|
745 "-*- texinfo -*-\n\ |
|
|
746 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ |
|
|
747 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ |
|
|
748 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ |
|
|
749 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ |
|
|
750 Return a matrix with Poisson distributed random elements. The arguments\n\ |
|
|
751 are handled the same as the arguments for @code{rand}, except for the\n\ |
|
|
752 argument @var{l}.\n\ |
|
|
753 \n\ |
|
|
754 Five different algorithms are used depending on the range of @var{l}\n\ |
|
|
755 and whether or not @var{l} is a scalar or a matrix.\n\ |
|
|
756 \n\ |
|
|
757 @table @asis\n\ |
|
|
758 @item For scalar @var{l} <= 12, use direct method.\n\ |
|
|
759 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
760 @item For scalar @var{l} > 12, use rejection method.[1]\n\ |
|
|
761 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
762 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ |
|
|
763 Stadlober E., et al., WinRand source code, available via FTP.\n\ |
|
|
764 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ |
|
|
765 Stadlober E., et al., WinRand source code, available via FTP, or\n\ |
|
|
766 H. Zechner, 'Efficient sampling from continuous and discrete\n\ |
|
|
767 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ |
|
|
768 University Graz, Austria, 1994.\n\ |
|
|
769 @item For @var{l} > 1e8, use normal approximation.\n\ |
|
|
770 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ |
|
|
771 D 50 p1284, 1994\n\ |
|
|
772 @end table\n\ |
|
|
773 @seealso{rand,randn,rande,randg}\n\ |
|
|
774 @end deftypefn") |
|
|
775 { |
|
|
776 octave_value retval; |
|
|
777 |
|
|
778 int nargin = args.length (); |
|
|
779 |
|
|
780 if (nargin < 1) |
|
|
781 error ("randp: insufficient arguments"); |
|
|
782 else |
|
|
783 { |
|
|
784 unwind_protect::begin_frame ("randp"); |
|
|
785 |
|
|
786 // This relies on the fact that elements are popped from the unwind |
|
|
787 // stack in the reverse of the order they are pushed |
|
|
788 // (i.e. current_distribution will be reset before calling |
|
|
789 // reset_rand_generator()). |
|
|
790 |
|
|
791 unwind_protect::add (reset_rand_generator, 0); |
|
|
792 unwind_protect_str (current_distribution); |
|
|
793 |
|
|
794 current_distribution = "poisson"; |
|
|
795 |
|
|
796 octave_rand::distribution (current_distribution); |
|
|
797 |
|
|
798 retval = do_rand (args, nargin, "randp", true); |
|
|
799 |
|
|
800 unwind_protect::run_frame ("randp"); |
|
|
801 } |
|
|
802 |
|
|
803 return retval; |
|
|
804 } |
|
|
805 |
|
|
806 /* |
|
|
807 %!test |
|
|
808 %! rand("state",12) |
|
|
809 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
810 %!test |
|
|
811 %! % statistical tests may fail occasionally. |
|
|
812 %! for a=[5 15] |
|
|
813 %! x = randp(a,100000,1); |
|
|
814 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
815 %! assert(mean(x),a,0.03); |
|
|
816 %! assert(var(x),a,0.2); |
|
|
817 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
818 %! assert(kurtosis(x),1/a,0.08); |
|
|
819 %! end |
|
|
820 %!test |
|
|
821 %! % statistical tests may fail occasionally. |
|
|
822 %! for a=[5 15] |
|
|
823 %! x = randp(a*ones(100000,1),100000,1); |
|
|
824 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
825 %! assert(mean(x),a,0.03); |
|
|
826 %! assert(var(x),a,0.2); |
|
|
827 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
828 %! assert(kurtosis(x),1/a,0.08); |
|
|
829 %! end |
|
|
830 %!test |
|
|
831 %! rand("seed",12) |
|
|
832 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
833 %!test |
|
|
834 %! % statistical tests may fail occasionally. |
|
|
835 %! for a=[5 15] |
|
|
836 %! x = randp(a,100000,1); |
|
|
837 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
838 %! assert(mean(x),a,0.03); |
|
|
839 %! assert(var(x),a,0.2); |
|
|
840 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
841 %! assert(kurtosis(x),1/a,0.08); |
|
|
842 %! end |
|
|
843 %!test |
|
|
844 %! % statistical tests may fail occasionally. |
|
|
845 %! for a=[5 15] |
|
|
846 %! x = randp(a*ones(100000,1),100000,1); |
|
|
847 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
848 %! assert(mean(x),a,0.03); |
|
|
849 %! assert(var(x),a,0.2); |
|
|
850 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
851 %! assert(kurtosis(x),1/a,0.08); |
|
|
852 %! end |
|
|
853 */ |
|
|
854 |
|
|
855 /* |
|
2928
|
856 ;;; Local Variables: *** |
|
|
857 ;;; mode: C++ *** |
|
|
858 ;;; End: *** |
|
|
859 */ |
