|
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 { |
|
5730
|
231 double d = args(idx+1).double_value (); |
|
2928
|
232 |
|
4543
|
233 if (! error_state) |
|
|
234 octave_rand::seed (d); |
|
|
235 } |
|
5730
|
236 else if (ts == "state" || ts == "twister") |
|
|
237 { |
|
|
238 ColumnVector s = |
|
|
239 ColumnVector (args(idx+1).vector_value(false, true)); |
|
|
240 |
|
|
241 if (! error_state) |
|
|
242 octave_rand::state (s); |
|
|
243 } |
|
4543
|
244 else |
|
4664
|
245 error ("%s: unrecognized string argument", fcn); |
|
4543
|
246 } |
|
|
247 else |
|
|
248 { |
|
|
249 dims.resize (nargin); |
|
|
250 |
|
|
251 for (int i = 0; i < nargin; i++) |
|
|
252 { |
|
5730
|
253 dims(i) = (octave_idx_type)args(idx+i).int_value (); |
|
4543
|
254 |
|
|
255 if (error_state) |
|
|
256 { |
|
4664
|
257 error ("%s: expecting integer arguments", fcn); |
|
4543
|
258 goto done; |
|
|
259 } |
|
|
260 } |
|
|
261 |
|
|
262 goto gen_matrix; |
|
|
263 } |
|
|
264 } |
|
|
265 break; |
|
2928
|
266 } |
|
|
267 |
|
4543
|
268 done: |
|
2928
|
269 |
|
|
270 return retval; |
|
|
271 |
|
|
272 gen_matrix: |
|
|
273 |
|
5355
|
274 dims.chop_trailing_singletons (); |
|
|
275 |
|
5730
|
276 if (additional_arg) |
|
|
277 { |
|
|
278 if (a.length() == 1) |
|
|
279 return octave_rand::nd_array (dims, a(0)); |
|
|
280 else |
|
|
281 { |
|
|
282 if (a.dims() != dims) |
|
|
283 { |
|
|
284 error ("%s: mismatch in argument size", fcn); |
|
|
285 return retval; |
|
|
286 } |
|
|
287 octave_idx_type len = a.length (); |
|
|
288 NDArray m (dims); |
|
|
289 double *v = m.fortran_vec (); |
|
|
290 for (octave_idx_type i = 0; i < len; i++) |
|
|
291 v[i] = octave_rand::scalar (a(i)); |
|
|
292 return m; |
|
|
293 } |
|
|
294 } |
|
|
295 else |
|
|
296 return octave_rand::nd_array (dims); |
|
2928
|
297 } |
|
|
298 |
|
4665
|
299 DEFUN_DLD (rand, args, , |
|
3369
|
300 "-*- texinfo -*-\n\ |
|
|
301 @deftypefn {Loadable Function} {} rand (@var{x})\n\ |
|
|
302 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ |
|
5730
|
303 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ |
|
|
304 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ |
|
3369
|
305 Return a matrix with random elements uniformly distributed on the\n\ |
|
|
306 interval (0, 1). The arguments are handled the same as the arguments\n\ |
|
5730
|
307 for @code{eye}.\n\ |
|
|
308 \n\ |
|
|
309 You can query the state of the random number generator using the\n\ |
|
3369
|
310 form\n\ |
|
2928
|
311 \n\ |
|
3369
|
312 @example\n\ |
|
5730
|
313 v = rand (\"state\")\n\ |
|
|
314 @end example\n\ |
|
|
315 \n\ |
|
|
316 This returns a column vector @var{v} of length 625. Later, you can\n\ |
|
|
317 restore the random number generator to the state @var{v}\n\ |
|
|
318 using the form\n\ |
|
|
319 \n\ |
|
|
320 @example\n\ |
|
|
321 rand (\"state\", v)\n\ |
|
3369
|
322 @end example\n\ |
|
|
323 \n\ |
|
|
324 @noindent\n\ |
|
5730
|
325 You may also initialize the state vector from an arbitrary vector of\n\ |
|
|
326 length <= 625 for @var{v}. This new state will be a hash based on the\n\ |
|
|
327 the value of @var{v}, not @var{v} itself.\n\ |
|
|
328 \n\ |
|
|
329 By default, the generator is initialized from @code{/dev/urandom} if it is\n\ |
|
|
330 available, otherwise from cpu time, wall clock time and the current\n\ |
|
|
331 fraction of a second.\n\ |
|
|
332 \n\ |
|
|
333 @code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\ |
|
|
334 (See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ |
|
|
335 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ |
|
|
336 Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\ |
|
|
337 @url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\ |
|
|
338 Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\ |
|
|
339 values together, otherwise the generator state can be learned after\n\ |
|
|
340 reading 624 consecutive values.\n\ |
|
|
341 \n\ |
|
|
342 @code{rand} includes a second random number generator, that was the\n\ |
|
|
343 previous generator used in octave. The new generator is used by default\n\ |
|
|
344 as it is significantly faster than the old generator, and produces\n\ |
|
|
345 random numebrs with a significantly longer cycle time. However, in\n\ |
|
|
346 some circumstances it might be desireable to obtain the same random\n\ |
|
|
347 sequences as used by the old generators. To do this the keyword\n\ |
|
|
348 \"seed\" is used to specify that the old generators should be use,\n\ |
|
|
349 as in\n\ |
|
2928
|
350 \n\ |
|
3369
|
351 @example\n\ |
|
5730
|
352 rand (\"seed\", val)\n\ |
|
3369
|
353 @end example\n\ |
|
|
354 \n\ |
|
5730
|
355 which sets the seed of the generator to @var{val}. The seed of the\n\ |
|
|
356 generator can be queried with\n\ |
|
|
357 \n\ |
|
|
358 @example\n\ |
|
|
359 s = rand (\"seed\")\n\ |
|
|
360 @end example\n\ |
|
|
361 \n\ |
|
|
362 However, it should be noted that querying the seed will not cause\n\ |
|
|
363 @code{rand} to use the old generators, only setting the seed will.\n\ |
|
|
364 To cause @code{rand} to once again use the new generators, the\n\ |
|
|
365 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ |
|
|
366 @seealso{randn,rande,randg,randp}\n\ |
|
3369
|
367 @end deftypefn") |
|
2928
|
368 { |
|
4307
|
369 octave_value retval; |
|
2928
|
370 |
|
|
371 int nargin = args.length (); |
|
|
372 |
|
4543
|
373 retval = do_rand (args, nargin, "rand"); |
|
2928
|
374 |
|
|
375 return retval; |
|
|
376 } |
|
|
377 |
|
5730
|
378 /* |
|
|
379 %!test # 'state' can be a scalar |
|
|
380 %! rand('state',12); x = rand(1,4); |
|
|
381 %! rand('state',12); y = rand(1,4); |
|
|
382 %! assert(x,y); |
|
|
383 %!test # 'state' can be a vector |
|
|
384 %! rand('state',[12,13]); x=rand(1,4); |
|
|
385 %! rand('state',[12;13]); y=rand(1,4); |
|
|
386 %! assert(x,y); |
|
|
387 %!test # querying 'state' doesn't disturb sequence |
|
|
388 %! rand('state',12); rand(1,2); x=rand(1,2); |
|
|
389 %! rand('state',12); rand(1,2); |
|
|
390 %! s=rand('state'); y=rand(1,2); |
|
|
391 %! assert(x,y); |
|
|
392 %! rand('state',s); z=rand(1,2); |
|
|
393 %! assert(x,z); |
|
|
394 %!test # 'seed' must be a scalar |
|
|
395 %! rand('seed',12); x = rand(1,4); |
|
|
396 %! rand('seed',12); y = rand(1,4); |
|
|
397 %! assert(x,y); |
|
|
398 %!error(rand('seed',[12,13])) |
|
|
399 %!test # querying 'seed' returns a value which can be used later |
|
|
400 %! s=rand('seed'); x=rand(1,2); |
|
|
401 %! rand('seed',s); y=rand(1,2); |
|
|
402 %! assert(x,y); |
|
|
403 %!test # querying 'seed' doesn't disturb sequence |
|
|
404 %! rand('seed',12); rand(1,2); x=rand(1,2); |
|
|
405 %! rand('seed',12); rand(1,2); |
|
|
406 %! s=rand('seed'); y=rand(1,2); |
|
|
407 %! assert(x,y); |
|
|
408 %! rand('seed',s); z=rand(1,2); |
|
|
409 %! assert(x,z); |
|
|
410 */ |
|
|
411 |
|
|
412 /* |
|
|
413 %!test |
|
|
414 %! % statistical tests may fail occasionally. |
|
|
415 %! rand("state",12); |
|
|
416 %! x = rand(100000,1); |
|
|
417 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
418 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
419 %! assert(mean(x),0.5,0.0024); |
|
|
420 %! assert(var(x),1/48,0.0632); |
|
|
421 %! assert(skewness(x),0,0.012); |
|
|
422 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
423 %!test |
|
|
424 %! % statistical tests may fail occasionally. |
|
|
425 %! rand("seed",12); |
|
|
426 %! x = rand(100000,1); |
|
|
427 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
428 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
429 %! assert(mean(x),0.5,0.0024); |
|
|
430 %! assert(var(x),1/48,0.0632); |
|
|
431 %! assert(skewness(x),0,0.012); |
|
|
432 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
433 */ |
|
|
434 |
|
|
435 |
|
4307
|
436 static std::string current_distribution = octave_rand::distribution (); |
|
|
437 |
|
2928
|
438 static void |
|
|
439 reset_rand_generator (void *) |
|
|
440 { |
|
4307
|
441 octave_rand::distribution (current_distribution); |
|
2928
|
442 } |
|
|
443 |
|
4665
|
444 DEFUN_DLD (randn, args, , |
|
3369
|
445 "-*- texinfo -*-\n\ |
|
|
446 @deftypefn {Loadable Function} {} randn (@var{x})\n\ |
|
|
447 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ |
|
5730
|
448 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
|
|
449 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ |
|
|
450 Return a matrix with normally distributed random elements. The\n\ |
|
|
451 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
3369
|
452 \n\ |
|
5730
|
453 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
454 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ |
|
|
455 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
456 J. Statistical Software, vol 5, 2000,\n\ |
|
|
457 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
2928
|
458 \n\ |
|
5730
|
459 @seealso{rand,rande,randg,randp}\n\ |
|
3369
|
460 @end deftypefn") |
|
2928
|
461 { |
|
4307
|
462 octave_value retval; |
|
2928
|
463 |
|
|
464 int nargin = args.length (); |
|
|
465 |
|
4543
|
466 unwind_protect::begin_frame ("randn"); |
|
2928
|
467 |
|
4543
|
468 // This relies on the fact that elements are popped from the unwind |
|
|
469 // stack in the reverse of the order they are pushed |
|
|
470 // (i.e. current_distribution will be reset before calling |
|
|
471 // reset_rand_generator()). |
|
2928
|
472 |
|
4543
|
473 unwind_protect::add (reset_rand_generator, 0); |
|
|
474 unwind_protect_str (current_distribution); |
|
2928
|
475 |
|
4543
|
476 current_distribution = "normal"; |
|
2928
|
477 |
|
4543
|
478 octave_rand::distribution (current_distribution); |
|
2928
|
479 |
|
4543
|
480 retval = do_rand (args, nargin, "randn"); |
|
2928
|
481 |
|
4543
|
482 unwind_protect::run_frame ("randn"); |
|
2928
|
483 |
|
|
484 return retval; |
|
|
485 } |
|
|
486 |
|
|
487 /* |
|
5730
|
488 %!test |
|
|
489 %! % statistical tests may fail occasionally. |
|
|
490 %! rand("state",12); |
|
|
491 %! x = randn(100000,1); |
|
|
492 %! assert(mean(x),0,0.01); |
|
|
493 %! assert(var(x),1,0.02); |
|
|
494 %! assert(skewness(x),0,0.02); |
|
|
495 %! assert(kurtosis(x),0,0.04); |
|
|
496 %!test |
|
|
497 %! % statistical tests may fail occasionally. |
|
|
498 %! rand("seed",12); |
|
|
499 %! x = randn(100000,1); |
|
|
500 %! assert(mean(x),0,0.01); |
|
|
501 %! assert(var(x),1,0.02); |
|
|
502 %! assert(skewness(x),0,0.02); |
|
|
503 %! assert(kurtosis(x),0,0.04); |
|
|
504 */ |
|
|
505 |
|
|
506 DEFUN_DLD (rande, args, , |
|
|
507 "-*- texinfo -*-\n\ |
|
|
508 @deftypefn {Loadable Function} {} rande (@var{x})\n\ |
|
|
509 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ |
|
|
510 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ |
|
|
511 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ |
|
|
512 Return a matrix with exponentially distributed random elements. The\n\ |
|
|
513 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
|
514 \n\ |
|
|
515 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
516 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ |
|
|
517 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
518 J. Statistical Software, vol 5, 2000,\n\ |
|
|
519 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
|
520 @seealso{rand,randn,randg,randp}\n\ |
|
|
521 @end deftypefn") |
|
|
522 { |
|
|
523 octave_value retval; |
|
|
524 |
|
|
525 int nargin = args.length (); |
|
|
526 |
|
|
527 unwind_protect::begin_frame ("rande"); |
|
|
528 |
|
|
529 // This relies on the fact that elements are popped from the unwind |
|
|
530 // stack in the reverse of the order they are pushed |
|
|
531 // (i.e. current_distribution will be reset before calling |
|
|
532 // reset_rand_generator()). |
|
|
533 |
|
|
534 unwind_protect::add (reset_rand_generator, 0); |
|
|
535 unwind_protect_str (current_distribution); |
|
|
536 |
|
|
537 current_distribution = "exponential"; |
|
|
538 |
|
|
539 octave_rand::distribution (current_distribution); |
|
|
540 |
|
|
541 retval = do_rand (args, nargin, "rande"); |
|
|
542 |
|
|
543 unwind_protect::run_frame ("rande"); |
|
|
544 |
|
|
545 return retval; |
|
|
546 } |
|
|
547 |
|
|
548 /* |
|
|
549 %!test |
|
|
550 %! % statistical tests may fail occasionally |
|
|
551 %! rand("state",12); |
|
|
552 %! x = rande(100000,1); |
|
|
553 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
554 %! assert(mean(x),1,0.01); |
|
|
555 %! assert(var(x),1,0.03); |
|
|
556 %! assert(skewness(x),2,0.06); |
|
|
557 %! assert(kurtosis(x),6,0.7); |
|
|
558 %!test |
|
|
559 %! % statistical tests may fail occasionally |
|
|
560 %! rand("seed",12); |
|
|
561 %! x = rande(100000,1); |
|
|
562 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
563 %! assert(mean(x),1,0.01); |
|
|
564 %! assert(var(x),1,0.03); |
|
|
565 %! assert(skewness(x),2,0.06); |
|
|
566 %! assert(kurtosis(x),6,0.7); |
|
|
567 */ |
|
|
568 |
|
|
569 DEFUN_DLD (randg, args, , |
|
|
570 "-*- texinfo -*-\n\ |
|
|
571 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ |
|
|
572 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ |
|
|
573 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ |
|
|
574 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ |
|
|
575 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ |
|
|
576 The arguments are handled the same as the arguments for @code{rand},\n\ |
|
|
577 except for the argument @var{a}.\n\ |
|
|
578 \n\ |
|
|
579 This can be used to generate many distributions:\n\ |
|
|
580 \n\ |
|
|
581 @table @asis\n\ |
|
|
582 @item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\ |
|
|
583 @example\n\ |
|
|
584 r = b*randg(a)\n\ |
|
|
585 @end example\n\ |
|
|
586 @item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\ |
|
|
587 @example\n\ |
|
|
588 r1 = randg(a,1)\n\ |
|
|
589 r = r1 / (r1 + randg(b,1))\n\ |
|
|
590 @end example\n\ |
|
|
591 @item @code{Erlang(a, n)}\n\ |
|
|
592 @example\n\ |
|
|
593 r = a*randg(n)\n\ |
|
|
594 @end example\n\ |
|
|
595 @item @code{chisq(df)} for @code{df > 0}\n\ |
|
|
596 @example\n\ |
|
|
597 r = 2*randg(df/2)\n\ |
|
|
598 @end example\n\ |
|
|
599 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ |
|
|
600 @example\n\ |
|
|
601 r = randn() / sqrt(2*randg(df/2)/df)\n\ |
|
|
602 @end example\n\ |
|
|
603 @item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\ |
|
|
604 @example\n\ |
|
|
605 r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\ |
|
|
606 r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\ |
|
|
607 r = r1 / r2\n\n\ |
|
|
608 @end example\n\ |
|
|
609 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ |
|
|
610 @example\n\ |
|
|
611 r = randp((1-p)/p * randg(n))\n\ |
|
|
612 @end example\n\ |
|
|
613 @item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\ |
|
|
614 (use chisq if @code{L = 0})\n\ |
|
|
615 @example\n\ |
|
|
616 r = randp(L/2)\n\ |
|
|
617 r(r > 0) = 2*randg(r(r > 0))\n\ |
|
|
618 r(df > 0) += 2*randg(df(df > 0)/2)\n\ |
|
|
619 @end example\n\ |
|
|
620 @item @code{Dirichlet(a1,...,ak)}\n\ |
|
|
621 @example\n\ |
|
|
622 r = (randg(a1),...,randg(ak))\n\ |
|
|
623 r = r / sum(r)\n\ |
|
|
624 @end example\n\ |
|
|
625 @end table\n\ |
|
|
626 @seealso{rand,randn,rande,randp}\n\ |
|
|
627 @end deftypefn") |
|
|
628 { |
|
|
629 octave_value retval; |
|
|
630 |
|
|
631 int nargin = args.length (); |
|
|
632 |
|
|
633 if (nargin < 1) |
|
|
634 error ("randg: insufficient arguments"); |
|
|
635 else |
|
|
636 { |
|
|
637 unwind_protect::begin_frame ("randg"); |
|
|
638 |
|
|
639 // This relies on the fact that elements are popped from the unwind |
|
|
640 // stack in the reverse of the order they are pushed |
|
|
641 // (i.e. current_distribution will be reset before calling |
|
|
642 // reset_rand_generator()). |
|
|
643 |
|
|
644 unwind_protect::add (reset_rand_generator, 0); |
|
|
645 unwind_protect_str (current_distribution); |
|
|
646 |
|
|
647 current_distribution = "gamma"; |
|
|
648 |
|
|
649 octave_rand::distribution (current_distribution); |
|
|
650 |
|
|
651 retval = do_rand (args, nargin, "randg", true); |
|
|
652 |
|
|
653 unwind_protect::run_frame ("randg"); |
|
|
654 } |
|
|
655 |
|
|
656 return retval; |
|
|
657 } |
|
|
658 |
|
|
659 /* |
|
|
660 %!test |
|
|
661 %! rand("state",12) |
|
|
662 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
663 %!test |
|
|
664 %! % statistical tests may fail occasionally. |
|
|
665 %! a=0.1; x = randg(a,100000,1); |
|
|
666 %! assert(mean(x), a, 0.01); |
|
|
667 %! assert(var(x), a, 0.01); |
|
|
668 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
669 %! assert(kurtosis(x),6/a, 50.); |
|
|
670 %!test |
|
|
671 %! % statistical tests may fail occasionally. |
|
|
672 %! a=0.95; x = randg(a,100000,1); |
|
|
673 %! assert(mean(x), a, 0.01); |
|
|
674 %! assert(var(x), a, 0.04); |
|
|
675 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
676 %! assert(kurtosis(x),6/a, 2.); |
|
|
677 %!test |
|
|
678 %! % statistical tests may fail occasionally. |
|
|
679 %! a=1; x = randg(a,100000,1); |
|
|
680 %! assert(mean(x),a, 0.01); |
|
|
681 %! assert(var(x),a, 0.04); |
|
|
682 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
683 %! assert(kurtosis(x),6/a, 2.); |
|
|
684 %!test |
|
|
685 %! % statistical tests may fail occasionally. |
|
|
686 %! a=10; x = randg(a,100000,1); |
|
|
687 %! assert(mean(x), a, 0.1); |
|
|
688 %! assert(var(x), a, 0.5); |
|
|
689 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
690 %! assert(kurtosis(x),6/a, 0.5); |
|
|
691 %!test |
|
|
692 %! % statistical tests may fail occasionally. |
|
|
693 %! a=100; x = randg(a,100000,1); |
|
|
694 %! assert(mean(x), a, 0.2); |
|
|
695 %! assert(var(x), a, 2.); |
|
|
696 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
697 %! assert(kurtosis(x),6/a, 0.2); |
|
|
698 %!test |
|
|
699 %! rand("seed",12) |
|
|
700 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
701 %!test |
|
|
702 %! % statistical tests may fail occasionally. |
|
|
703 %! a=0.1; x = randg(a,100000,1); |
|
|
704 %! assert(mean(x), a, 0.01); |
|
|
705 %! assert(var(x), a, 0.01); |
|
|
706 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
707 %! assert(kurtosis(x),6/a, 50.); |
|
|
708 %!test |
|
|
709 %! % statistical tests may fail occasionally. |
|
|
710 %! a=0.95; x = randg(a,100000,1); |
|
|
711 %! assert(mean(x), a, 0.01); |
|
|
712 %! assert(var(x), a, 0.04); |
|
|
713 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
714 %! assert(kurtosis(x),6/a, 2.); |
|
|
715 %!test |
|
|
716 %! % statistical tests may fail occasionally. |
|
|
717 %! a=1; x = randg(a,100000,1); |
|
|
718 %! assert(mean(x),a, 0.01); |
|
|
719 %! assert(var(x),a, 0.04); |
|
|
720 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
721 %! assert(kurtosis(x),6/a, 2.); |
|
|
722 %!test |
|
|
723 %! % statistical tests may fail occasionally. |
|
|
724 %! a=10; x = randg(a,100000,1); |
|
|
725 %! assert(mean(x), a, 0.1); |
|
|
726 %! assert(var(x), a, 0.5); |
|
|
727 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
728 %! assert(kurtosis(x),6/a, 0.5); |
|
|
729 %!test |
|
|
730 %! % statistical tests may fail occasionally. |
|
|
731 %! a=100; x = randg(a,100000,1); |
|
|
732 %! assert(mean(x), a, 0.2); |
|
|
733 %! assert(var(x), a, 2.); |
|
|
734 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
735 %! assert(kurtosis(x),6/a, 0.2); |
|
|
736 */ |
|
|
737 |
|
|
738 |
|
|
739 DEFUN_DLD (randp, args, , |
|
|
740 "-*- texinfo -*-\n\ |
|
|
741 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ |
|
|
742 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ |
|
|
743 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ |
|
|
744 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ |
|
|
745 Return a matrix with Poisson distributed random elements. The arguments\n\ |
|
|
746 are handled the same as the arguments for @code{rand}, except for the\n\ |
|
|
747 argument @var{l}.\n\ |
|
|
748 \n\ |
|
|
749 Five different algorithms are used depending on the range of @var{l}\n\ |
|
|
750 and whether or not @var{l} is a scalar or a matrix.\n\ |
|
|
751 \n\ |
|
|
752 @table @asis\n\ |
|
|
753 @item For scalar @var{l} <= 12, use direct method.\n\ |
|
|
754 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
755 @item For scalar @var{l} > 12, use rejection method.[1]\n\ |
|
|
756 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
757 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ |
|
|
758 Stadlober E., et al., WinRand source code, available via FTP.\n\ |
|
|
759 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ |
|
|
760 Stadlober E., et al., WinRand source code, available via FTP, or\n\ |
|
|
761 H. Zechner, 'Efficient sampling from continuous and discrete\n\ |
|
|
762 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ |
|
|
763 University Graz, Austria, 1994.\n\ |
|
|
764 @item For @var{l} > 1e8, use normal approximation.\n\ |
|
|
765 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ |
|
|
766 D 50 p1284, 1994\n\ |
|
|
767 @end table\n\ |
|
|
768 @seealso{rand,randn,rande,randg}\n\ |
|
|
769 @end deftypefn") |
|
|
770 { |
|
|
771 octave_value retval; |
|
|
772 |
|
|
773 int nargin = args.length (); |
|
|
774 |
|
|
775 if (nargin < 1) |
|
|
776 error ("randp: insufficient arguments"); |
|
|
777 else |
|
|
778 { |
|
|
779 unwind_protect::begin_frame ("randp"); |
|
|
780 |
|
|
781 // This relies on the fact that elements are popped from the unwind |
|
|
782 // stack in the reverse of the order they are pushed |
|
|
783 // (i.e. current_distribution will be reset before calling |
|
|
784 // reset_rand_generator()). |
|
|
785 |
|
|
786 unwind_protect::add (reset_rand_generator, 0); |
|
|
787 unwind_protect_str (current_distribution); |
|
|
788 |
|
|
789 current_distribution = "poisson"; |
|
|
790 |
|
|
791 octave_rand::distribution (current_distribution); |
|
|
792 |
|
|
793 retval = do_rand (args, nargin, "randp", true); |
|
|
794 |
|
|
795 unwind_protect::run_frame ("randp"); |
|
|
796 } |
|
|
797 |
|
|
798 return retval; |
|
|
799 } |
|
|
800 |
|
|
801 /* |
|
|
802 %!test |
|
|
803 %! rand("state",12) |
|
|
804 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
805 %!test |
|
|
806 %! % statistical tests may fail occasionally. |
|
|
807 %! for a=[5 15] |
|
|
808 %! x = randp(a,100000,1); |
|
|
809 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
810 %! assert(mean(x),a,0.03); |
|
|
811 %! assert(var(x),a,0.2); |
|
|
812 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
813 %! assert(kurtosis(x),1/a,0.08); |
|
|
814 %! end |
|
|
815 %!test |
|
|
816 %! % statistical tests may fail occasionally. |
|
|
817 %! for a=[5 15] |
|
|
818 %! x = randp(a*ones(100000,1),100000,1); |
|
|
819 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
820 %! assert(mean(x),a,0.03); |
|
|
821 %! assert(var(x),a,0.2); |
|
|
822 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
823 %! assert(kurtosis(x),1/a,0.08); |
|
|
824 %! end |
|
|
825 %!test |
|
|
826 %! rand("seed",12) |
|
|
827 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
828 %!test |
|
|
829 %! % statistical tests may fail occasionally. |
|
|
830 %! for a=[5 15] |
|
|
831 %! x = randp(a,100000,1); |
|
|
832 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
833 %! assert(mean(x),a,0.03); |
|
|
834 %! assert(var(x),a,0.2); |
|
|
835 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
836 %! assert(kurtosis(x),1/a,0.08); |
|
|
837 %! end |
|
|
838 %!test |
|
|
839 %! % statistical tests may fail occasionally. |
|
|
840 %! for a=[5 15] |
|
|
841 %! x = randp(a*ones(100000,1),100000,1); |
|
|
842 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
843 %! assert(mean(x),a,0.03); |
|
|
844 %! assert(var(x),a,0.2); |
|
|
845 %! assert(skewness(x),1/sqrt(a),0.03); |
|
|
846 %! assert(kurtosis(x),1/a,0.08); |
|
|
847 %! end |
|
|
848 */ |
|
|
849 |
|
|
850 /* |
|
2928
|
851 ;;; Local Variables: *** |
|
|
852 ;;; mode: C++ *** |
|
|
853 ;;; End: *** |
|
|
854 */ |
