|
2928
|
1 /* |
|
|
2 |
|
7017
|
3 Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2003, 2005, 2006, |
|
|
4 2007 John W. Eaton |
|
2928
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
7016
|
10 Free Software Foundation; either version 3 of the License, or (at your |
|
|
11 option) any later version. |
|
2928
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
7016
|
19 along with Octave; see the file COPYING. If not, see |
|
|
20 <http://www.gnu.org/licenses/>. |
|
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 |
|
6437
|
44 /* |
|
|
45 %!shared __random_statistical_tests__ |
|
|
46 %! % Flag whether the statistical tests should be run in "make check" or not |
|
|
47 %! __random_statistical_tests__ = 0; |
|
|
48 */ |
|
|
49 |
|
4307
|
50 static octave_value |
|
5730
|
51 do_rand (const octave_value_list& args, int nargin, const char *fcn, |
|
|
52 bool additional_arg = false) |
|
2928
|
53 { |
|
4307
|
54 octave_value retval; |
|
5730
|
55 NDArray a; |
|
|
56 int idx = 0; |
|
|
57 dim_vector dims; |
|
2928
|
58 |
|
5730
|
59 if (additional_arg) |
|
|
60 { |
|
|
61 if (nargin == 0) |
|
|
62 { |
|
|
63 error ("%s: expecting at least one argument", fcn); |
|
|
64 goto done; |
|
|
65 } |
|
|
66 else if (args(0).is_string()) |
|
|
67 additional_arg = false; |
|
|
68 else |
|
|
69 { |
|
|
70 a = args(0).array_value (); |
|
|
71 if (error_state) |
|
|
72 { |
|
|
73 error ("%s: expecting scalar or matrix arguments", fcn); |
|
|
74 goto done; |
|
|
75 } |
|
|
76 idx++; |
|
|
77 nargin--; |
|
|
78 } |
|
|
79 } |
|
2928
|
80 |
|
4543
|
81 switch (nargin) |
|
2928
|
82 { |
|
4543
|
83 case 0: |
|
|
84 { |
|
5730
|
85 if (additional_arg) |
|
|
86 dims = a.dims (); |
|
|
87 else |
|
|
88 { |
|
|
89 dims.resize (2); |
|
4543
|
90 |
|
5730
|
91 dims(0) = 1; |
|
|
92 dims(1) = 1; |
|
|
93 } |
|
4543
|
94 goto gen_matrix; |
|
|
95 } |
|
|
96 break; |
|
2928
|
97 |
|
4543
|
98 case 1: |
|
|
99 { |
|
5730
|
100 octave_value tmp = args(idx); |
|
4543
|
101 |
|
|
102 if (tmp.is_string ()) |
|
|
103 { |
|
|
104 std::string s_arg = tmp.string_value (); |
|
2928
|
105 |
|
4543
|
106 if (s_arg == "dist") |
|
|
107 { |
|
|
108 retval = octave_rand::distribution (); |
|
|
109 } |
|
5730
|
110 else if (s_arg == "seed") |
|
4543
|
111 { |
|
|
112 retval = octave_rand::seed (); |
|
|
113 } |
|
5730
|
114 else if (s_arg == "state" || s_arg == "twister") |
|
|
115 { |
|
|
116 retval = octave_rand::state (); |
|
|
117 } |
|
4543
|
118 else if (s_arg == "uniform") |
|
|
119 { |
|
|
120 octave_rand::uniform_distribution (); |
|
|
121 } |
|
|
122 else if (s_arg == "normal") |
|
|
123 { |
|
|
124 octave_rand::normal_distribution (); |
|
|
125 } |
|
5730
|
126 else if (s_arg == "exponential") |
|
|
127 { |
|
|
128 octave_rand::exponential_distribution (); |
|
|
129 } |
|
|
130 else if (s_arg == "poisson") |
|
|
131 { |
|
|
132 octave_rand::poisson_distribution (); |
|
|
133 } |
|
|
134 else if (s_arg == "gamma") |
|
|
135 { |
|
|
136 octave_rand::gamma_distribution (); |
|
|
137 } |
|
4543
|
138 else |
|
4664
|
139 error ("%s: unrecognized string argument", fcn); |
|
4543
|
140 } |
|
|
141 else if (tmp.is_scalar_type ()) |
|
|
142 { |
|
|
143 double dval = tmp.double_value (); |
|
2928
|
144 |
|
4543
|
145 if (xisnan (dval)) |
|
|
146 { |
|
4664
|
147 error ("%s: NaN is invalid a matrix dimension", fcn); |
|
4543
|
148 } |
|
|
149 else |
|
|
150 { |
|
|
151 dims.resize (2); |
|
|
152 |
|
5275
|
153 dims(0) = NINTbig (tmp.double_value ()); |
|
|
154 dims(1) = NINTbig (tmp.double_value ()); |
|
2928
|
155 |
|
4543
|
156 if (! error_state) |
|
|
157 goto gen_matrix; |
|
|
158 } |
|
|
159 } |
|
|
160 else if (tmp.is_range ()) |
|
|
161 { |
|
|
162 Range r = tmp.range_value (); |
|
|
163 |
|
|
164 if (r.all_elements_are_ints ()) |
|
|
165 { |
|
5275
|
166 octave_idx_type n = r.nelem (); |
|
4543
|
167 |
|
|
168 dims.resize (n); |
|
|
169 |
|
5275
|
170 octave_idx_type base = NINTbig (r.base ()); |
|
|
171 octave_idx_type incr = NINTbig (r.inc ()); |
|
|
172 octave_idx_type lim = NINTbig (r.limit ()); |
|
2928
|
173 |
|
4543
|
174 if (base < 0 || lim < 0) |
|
4664
|
175 error ("%s: all dimensions must be nonnegative", fcn); |
|
4543
|
176 else |
|
|
177 { |
|
5275
|
178 for (octave_idx_type i = 0; i < n; i++) |
|
4543
|
179 { |
|
|
180 dims(i) = base; |
|
|
181 base += incr; |
|
|
182 } |
|
2928
|
183 |
|
4543
|
184 goto gen_matrix; |
|
|
185 } |
|
|
186 } |
|
|
187 else |
|
4664
|
188 error ("%s: expecting all elements of range to be integers", |
|
|
189 fcn); |
|
4543
|
190 } |
|
|
191 else if (tmp.is_matrix_type ()) |
|
|
192 { |
|
|
193 Array<int> iv = tmp.int_vector_value (true); |
|
|
194 |
|
|
195 if (! error_state) |
|
|
196 { |
|
5275
|
197 octave_idx_type len = iv.length (); |
|
2928
|
198 |
|
4543
|
199 dims.resize (len); |
|
|
200 |
|
5275
|
201 for (octave_idx_type i = 0; i < len; i++) |
|
4543
|
202 { |
|
5275
|
203 octave_idx_type elt = iv(i); |
|
4543
|
204 |
|
|
205 if (elt < 0) |
|
|
206 { |
|
4664
|
207 error ("%s: all dimensions must be nonnegative", fcn); |
|
4543
|
208 goto done; |
|
|
209 } |
|
|
210 |
|
|
211 dims(i) = iv(i); |
|
|
212 } |
|
2928
|
213 |
|
4543
|
214 goto gen_matrix; |
|
|
215 } |
|
|
216 else |
|
4664
|
217 error ("%s: expecting integer vector", fcn); |
|
4543
|
218 } |
|
|
219 else |
|
|
220 { |
|
|
221 gripe_wrong_type_arg ("rand", tmp); |
|
|
222 return retval; |
|
|
223 } |
|
|
224 } |
|
|
225 break; |
|
|
226 |
|
|
227 default: |
|
|
228 { |
|
5730
|
229 octave_value tmp = args(idx); |
|
4543
|
230 |
|
|
231 if (nargin == 2 && tmp.is_string ()) |
|
|
232 { |
|
5164
|
233 std::string ts = tmp.string_value (); |
|
|
234 |
|
5730
|
235 if (ts == "seed") |
|
4543
|
236 { |
|
5782
|
237 if (args(idx+1).is_real_scalar ()) |
|
|
238 { |
|
|
239 double d = args(idx+1).double_value (); |
|
2928
|
240 |
|
5782
|
241 if (! error_state) |
|
|
242 octave_rand::seed (d); |
|
|
243 } |
|
|
244 else |
|
|
245 error ("%s: seed must be a real scalar", fcn); |
|
4543
|
246 } |
|
5730
|
247 else if (ts == "state" || ts == "twister") |
|
|
248 { |
|
|
249 ColumnVector s = |
|
|
250 ColumnVector (args(idx+1).vector_value(false, true)); |
|
|
251 |
|
|
252 if (! error_state) |
|
|
253 octave_rand::state (s); |
|
|
254 } |
|
4543
|
255 else |
|
4664
|
256 error ("%s: unrecognized string argument", fcn); |
|
4543
|
257 } |
|
|
258 else |
|
|
259 { |
|
|
260 dims.resize (nargin); |
|
|
261 |
|
|
262 for (int i = 0; i < nargin; i++) |
|
|
263 { |
|
5760
|
264 dims(i) = args(idx+i).int_value (); |
|
4543
|
265 |
|
|
266 if (error_state) |
|
|
267 { |
|
4664
|
268 error ("%s: expecting integer arguments", fcn); |
|
4543
|
269 goto done; |
|
|
270 } |
|
|
271 } |
|
|
272 |
|
|
273 goto gen_matrix; |
|
|
274 } |
|
|
275 } |
|
|
276 break; |
|
2928
|
277 } |
|
|
278 |
|
4543
|
279 done: |
|
2928
|
280 |
|
|
281 return retval; |
|
|
282 |
|
|
283 gen_matrix: |
|
|
284 |
|
5355
|
285 dims.chop_trailing_singletons (); |
|
|
286 |
|
5730
|
287 if (additional_arg) |
|
|
288 { |
|
|
289 if (a.length() == 1) |
|
|
290 return octave_rand::nd_array (dims, a(0)); |
|
|
291 else |
|
|
292 { |
|
|
293 if (a.dims() != dims) |
|
|
294 { |
|
|
295 error ("%s: mismatch in argument size", fcn); |
|
|
296 return retval; |
|
|
297 } |
|
|
298 octave_idx_type len = a.length (); |
|
|
299 NDArray m (dims); |
|
|
300 double *v = m.fortran_vec (); |
|
|
301 for (octave_idx_type i = 0; i < len; i++) |
|
|
302 v[i] = octave_rand::scalar (a(i)); |
|
|
303 return m; |
|
|
304 } |
|
|
305 } |
|
|
306 else |
|
|
307 return octave_rand::nd_array (dims); |
|
2928
|
308 } |
|
|
309 |
|
4665
|
310 DEFUN_DLD (rand, args, , |
|
3369
|
311 "-*- texinfo -*-\n\ |
|
|
312 @deftypefn {Loadable Function} {} rand (@var{x})\n\ |
|
|
313 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\ |
|
5730
|
314 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\ |
|
|
315 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\ |
|
3369
|
316 Return a matrix with random elements uniformly distributed on the\n\ |
|
|
317 interval (0, 1). The arguments are handled the same as the arguments\n\ |
|
5730
|
318 for @code{eye}.\n\ |
|
|
319 \n\ |
|
|
320 You can query the state of the random number generator using the\n\ |
|
3369
|
321 form\n\ |
|
2928
|
322 \n\ |
|
3369
|
323 @example\n\ |
|
5730
|
324 v = rand (\"state\")\n\ |
|
|
325 @end example\n\ |
|
|
326 \n\ |
|
|
327 This returns a column vector @var{v} of length 625. Later, you can\n\ |
|
|
328 restore the random number generator to the state @var{v}\n\ |
|
|
329 using the form\n\ |
|
|
330 \n\ |
|
|
331 @example\n\ |
|
|
332 rand (\"state\", v)\n\ |
|
3369
|
333 @end example\n\ |
|
|
334 \n\ |
|
|
335 @noindent\n\ |
|
5730
|
336 You may also initialize the state vector from an arbitrary vector of\n\ |
|
|
337 length <= 625 for @var{v}. This new state will be a hash based on the\n\ |
|
5798
|
338 value of @var{v}, not @var{v} itself.\n\ |
|
5730
|
339 \n\ |
|
|
340 By default, the generator is initialized from @code{/dev/urandom} if it is\n\ |
|
|
341 available, otherwise from cpu time, wall clock time and the current\n\ |
|
|
342 fraction of a second.\n\ |
|
|
343 \n\ |
|
7096
|
344 To compute the psuedo-random sequence, @code{rand} uses the Mersenne\n\ |
|
|
345 Twister with a period of 2^19937-1 (See M. Matsumoto and T. Nishimura,\n\ |
|
|
346 ``Mersenne Twister: A 623-dimensionally\n\ |
|
5730
|
347 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ |
|
7001
|
348 Modeling and Computer Simulation Vol. 8, No. 1, January pp.3-30 1998,\n\ |
|
7171
|
349 @url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}).\n\ |
|
6547
|
350 Do @strong{not} use for cryptography without securely hashing\n\ |
|
|
351 several returned values together, otherwise the generator state\n\ |
|
|
352 can be learned after reading 624 consecutive values.\n\ |
|
5730
|
353 \n\ |
|
7096
|
354 Older versions of Octave used a different random number generator.\n\ |
|
|
355 The new generator is used by default\n\ |
|
5730
|
356 as it is significantly faster than the old generator, and produces\n\ |
|
5798
|
357 random numbers with a significantly longer cycle time. However, in\n\ |
|
|
358 some circumstances it might be desirable to obtain the same random\n\ |
|
5730
|
359 sequences as used by the old generators. To do this the keyword\n\ |
|
|
360 \"seed\" is used to specify that the old generators should be use,\n\ |
|
|
361 as in\n\ |
|
2928
|
362 \n\ |
|
3369
|
363 @example\n\ |
|
5730
|
364 rand (\"seed\", val)\n\ |
|
3369
|
365 @end example\n\ |
|
|
366 \n\ |
|
5730
|
367 which sets the seed of the generator to @var{val}. The seed of the\n\ |
|
|
368 generator can be queried with\n\ |
|
|
369 \n\ |
|
|
370 @example\n\ |
|
|
371 s = rand (\"seed\")\n\ |
|
|
372 @end example\n\ |
|
|
373 \n\ |
|
|
374 However, it should be noted that querying the seed will not cause\n\ |
|
|
375 @code{rand} to use the old generators, only setting the seed will.\n\ |
|
|
376 To cause @code{rand} to once again use the new generators, the\n\ |
|
|
377 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ |
|
5798
|
378 @seealso{randn, rande, randg, randp}\n\ |
|
3369
|
379 @end deftypefn") |
|
2928
|
380 { |
|
4307
|
381 octave_value retval; |
|
2928
|
382 |
|
|
383 int nargin = args.length (); |
|
|
384 |
|
4543
|
385 retval = do_rand (args, nargin, "rand"); |
|
2928
|
386 |
|
|
387 return retval; |
|
|
388 } |
|
|
389 |
|
5730
|
390 /* |
|
|
391 %!test # 'state' can be a scalar |
|
|
392 %! rand('state',12); x = rand(1,4); |
|
|
393 %! rand('state',12); y = rand(1,4); |
|
|
394 %! assert(x,y); |
|
|
395 %!test # 'state' can be a vector |
|
|
396 %! rand('state',[12,13]); x=rand(1,4); |
|
|
397 %! rand('state',[12;13]); y=rand(1,4); |
|
|
398 %! assert(x,y); |
|
|
399 %!test # querying 'state' doesn't disturb sequence |
|
|
400 %! rand('state',12); rand(1,2); x=rand(1,2); |
|
|
401 %! rand('state',12); rand(1,2); |
|
|
402 %! s=rand('state'); y=rand(1,2); |
|
|
403 %! assert(x,y); |
|
|
404 %! rand('state',s); z=rand(1,2); |
|
|
405 %! assert(x,z); |
|
|
406 %!test # 'seed' must be a scalar |
|
|
407 %! rand('seed',12); x = rand(1,4); |
|
|
408 %! rand('seed',12); y = rand(1,4); |
|
|
409 %! assert(x,y); |
|
|
410 %!error(rand('seed',[12,13])) |
|
|
411 %!test # querying 'seed' returns a value which can be used later |
|
|
412 %! s=rand('seed'); x=rand(1,2); |
|
|
413 %! rand('seed',s); y=rand(1,2); |
|
|
414 %! assert(x,y); |
|
|
415 %!test # querying 'seed' doesn't disturb sequence |
|
|
416 %! rand('seed',12); rand(1,2); x=rand(1,2); |
|
|
417 %! rand('seed',12); rand(1,2); |
|
|
418 %! s=rand('seed'); y=rand(1,2); |
|
|
419 %! assert(x,y); |
|
|
420 %! rand('seed',s); z=rand(1,2); |
|
|
421 %! assert(x,z); |
|
|
422 */ |
|
|
423 |
|
|
424 /* |
|
|
425 %!test |
|
6437
|
426 %! % Test fixed state |
|
|
427 %! rand("state",1); |
|
6443
|
428 %! assert (rand(1,6), [0.1343642441124013 0.8474337369372327 0.763774618976614 0.2550690257394218 0.495435087091941 0.4494910647887382],1e-6); |
|
6437
|
429 %!test |
|
6443
|
430 %! % Test fixed seed |
|
6437
|
431 %! rand("seed",1); |
|
6443
|
432 %! assert (rand(1,6), [0.8668024251237512 0.9126510815694928 0.09366085007786751 0.1664607301354408 0.7408077004365623 0.7615650338120759],1e-6); |
|
5730
|
433 %!test |
|
6437
|
434 %! if (__random_statistical_tests__) |
|
|
435 %! % statistical tests may fail occasionally. |
|
|
436 %! rand("state",12); |
|
|
437 %! x = rand(100000,1); |
|
|
438 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
439 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
440 %! assert(mean(x),0.5,0.0024); |
|
|
441 %! assert(var(x),1/48,0.0632); |
|
|
442 %! assert(skewness(x),0,0.012); |
|
|
443 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
444 %! endif |
|
|
445 %!test |
|
|
446 %! if (__random_statistical_tests__) |
|
|
447 %! % statistical tests may fail occasionally. |
|
|
448 %! rand("seed",12); |
|
|
449 %! x = rand(100000,1); |
|
|
450 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
451 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
452 %! assert(mean(x),0.5,0.0024); |
|
|
453 %! assert(var(x),1/48,0.0632); |
|
|
454 %! assert(skewness(x),0,0.012); |
|
|
455 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
456 %! endif |
|
5730
|
457 */ |
|
|
458 |
|
|
459 |
|
4307
|
460 static std::string current_distribution = octave_rand::distribution (); |
|
|
461 |
|
2928
|
462 static void |
|
|
463 reset_rand_generator (void *) |
|
|
464 { |
|
4307
|
465 octave_rand::distribution (current_distribution); |
|
2928
|
466 } |
|
|
467 |
|
4665
|
468 DEFUN_DLD (randn, args, , |
|
3369
|
469 "-*- texinfo -*-\n\ |
|
|
470 @deftypefn {Loadable Function} {} randn (@var{x})\n\ |
|
|
471 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ |
|
5730
|
472 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
|
|
473 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ |
|
|
474 Return a matrix with normally distributed random elements. The\n\ |
|
|
475 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
3369
|
476 \n\ |
|
5730
|
477 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
478 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ |
|
|
479 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
480 J. Statistical Software, vol 5, 2000,\n\ |
|
|
481 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
2928
|
482 \n\ |
|
6547
|
483 @seealso{rand, rande, randg, randp}\n\ |
|
3369
|
484 @end deftypefn") |
|
2928
|
485 { |
|
4307
|
486 octave_value retval; |
|
2928
|
487 |
|
|
488 int nargin = args.length (); |
|
|
489 |
|
4543
|
490 unwind_protect::begin_frame ("randn"); |
|
2928
|
491 |
|
4543
|
492 // This relies on the fact that elements are popped from the unwind |
|
|
493 // stack in the reverse of the order they are pushed |
|
|
494 // (i.e. current_distribution will be reset before calling |
|
|
495 // reset_rand_generator()). |
|
2928
|
496 |
|
4543
|
497 unwind_protect::add (reset_rand_generator, 0); |
|
|
498 unwind_protect_str (current_distribution); |
|
2928
|
499 |
|
4543
|
500 current_distribution = "normal"; |
|
2928
|
501 |
|
4543
|
502 octave_rand::distribution (current_distribution); |
|
2928
|
503 |
|
4543
|
504 retval = do_rand (args, nargin, "randn"); |
|
2928
|
505 |
|
4543
|
506 unwind_protect::run_frame ("randn"); |
|
2928
|
507 |
|
|
508 return retval; |
|
|
509 } |
|
|
510 |
|
|
511 /* |
|
5730
|
512 %!test |
|
6437
|
513 %! % Test fixed state |
|
|
514 %! randn("state",1); |
|
6443
|
515 %! assert (randn(1,6), [-2.666521678978671 -0.7381719971724564 1.507903992673601 0.6019427189162239 -0.450661261143348 -0.7054431351574116],1e-6); |
|
6437
|
516 %!test |
|
6443
|
517 %! % Test fixed seed |
|
6437
|
518 %! randn("seed",1); |
|
6443
|
519 %! assert (randn(1,6), [-1.039402365684509 -1.25938892364502 0.1968704611063004 0.3874166905879974 -0.5976632833480835 -0.6615074276924133],1e-6); |
|
5730
|
520 %!test |
|
6437
|
521 %! if (__random_statistical_tests__) |
|
|
522 %! % statistical tests may fail occasionally. |
|
|
523 %! randn("state",12); |
|
|
524 %! x = randn(100000,1); |
|
|
525 %! assert(mean(x),0,0.01); |
|
|
526 %! assert(var(x),1,0.02); |
|
|
527 %! assert(skewness(x),0,0.02); |
|
|
528 %! assert(kurtosis(x),0,0.04); |
|
|
529 %! endif |
|
|
530 %!test |
|
|
531 %! if (__random_statistical_tests__) |
|
|
532 %! % statistical tests may fail occasionally. |
|
|
533 %! randn("seed",12); |
|
|
534 %! x = randn(100000,1); |
|
|
535 %! assert(mean(x),0,0.01); |
|
|
536 %! assert(var(x),1,0.02); |
|
|
537 %! assert(skewness(x),0,0.02); |
|
|
538 %! assert(kurtosis(x),0,0.04); |
|
|
539 %! endif |
|
5730
|
540 */ |
|
|
541 |
|
|
542 DEFUN_DLD (rande, args, , |
|
|
543 "-*- texinfo -*-\n\ |
|
|
544 @deftypefn {Loadable Function} {} rande (@var{x})\n\ |
|
|
545 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ |
|
|
546 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ |
|
|
547 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ |
|
|
548 Return a matrix with exponentially distributed random elements. The\n\ |
|
|
549 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
|
550 \n\ |
|
|
551 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
552 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ |
|
|
553 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
554 J. Statistical Software, vol 5, 2000,\n\ |
|
|
555 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
6547
|
556 @seealso{rand, randn, randg, randp}\n\ |
|
5730
|
557 @end deftypefn") |
|
|
558 { |
|
|
559 octave_value retval; |
|
|
560 |
|
|
561 int nargin = args.length (); |
|
|
562 |
|
|
563 unwind_protect::begin_frame ("rande"); |
|
|
564 |
|
|
565 // This relies on the fact that elements are popped from the unwind |
|
|
566 // stack in the reverse of the order they are pushed |
|
|
567 // (i.e. current_distribution will be reset before calling |
|
|
568 // reset_rand_generator()). |
|
|
569 |
|
|
570 unwind_protect::add (reset_rand_generator, 0); |
|
|
571 unwind_protect_str (current_distribution); |
|
|
572 |
|
|
573 current_distribution = "exponential"; |
|
|
574 |
|
|
575 octave_rand::distribution (current_distribution); |
|
|
576 |
|
|
577 retval = do_rand (args, nargin, "rande"); |
|
|
578 |
|
|
579 unwind_protect::run_frame ("rande"); |
|
|
580 |
|
|
581 return retval; |
|
|
582 } |
|
|
583 |
|
|
584 /* |
|
|
585 %!test |
|
6437
|
586 %! % Test fixed state |
|
|
587 %! rande("state",1); |
|
6443
|
588 %! assert (rande(1,6), [3.602973885835625 0.1386190677555021 0.6743112889616958 0.4512830847258422 0.7255744741233175 0.3415969205292291],1e-6); |
|
6437
|
589 %!test |
|
6443
|
590 %! % Test fixed seed |
|
6437
|
591 %! rande("seed",1); |
|
6443
|
592 %! assert (rande(1,6), [0.06492075175653866 1.717980206012726 0.4816154008731246 0.5231300676241517 0.103910739364359 1.668931916356087],1e-6); |
|
5730
|
593 %!test |
|
6437
|
594 %! if (__random_statistical_tests__) |
|
|
595 %! % statistical tests may fail occasionally |
|
|
596 %! rande("state",1); |
|
|
597 %! x = rande(100000,1); |
|
|
598 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
599 %! assert(mean(x),1,0.01); |
|
|
600 %! assert(var(x),1,0.03); |
|
|
601 %! assert(skewness(x),2,0.06); |
|
|
602 %! assert(kurtosis(x),6,0.7); |
|
|
603 %! endif |
|
|
604 %!test |
|
|
605 %! if (__random_statistical_tests__) |
|
|
606 %! % statistical tests may fail occasionally |
|
|
607 %! rande("seed",1); |
|
|
608 %! x = rande(100000,1); |
|
|
609 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
610 %! assert(mean(x),1,0.01); |
|
|
611 %! assert(var(x),1,0.03); |
|
|
612 %! assert(skewness(x),2,0.06); |
|
|
613 %! assert(kurtosis(x),6,0.7); |
|
|
614 %! endif |
|
5730
|
615 */ |
|
|
616 |
|
|
617 DEFUN_DLD (randg, args, , |
|
|
618 "-*- texinfo -*-\n\ |
|
|
619 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ |
|
|
620 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ |
|
|
621 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ |
|
|
622 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ |
|
|
623 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ |
|
|
624 The arguments are handled the same as the arguments for @code{rand},\n\ |
|
|
625 except for the argument @var{a}.\n\ |
|
|
626 \n\ |
|
|
627 This can be used to generate many distributions:\n\ |
|
|
628 \n\ |
|
|
629 @table @asis\n\ |
|
6547
|
630 @item @code{gamma (a, b)} for @code{a > -1}, @code{b > 0}\n\ |
|
5730
|
631 @example\n\ |
|
6547
|
632 r = b * randg (a)\n\ |
|
5730
|
633 @end example\n\ |
|
6547
|
634 @item @code{beta (a, b)} for @code{a > -1}, @code{b > -1}\n\ |
|
5730
|
635 @example\n\ |
|
6547
|
636 r1 = randg (a, 1)\n\ |
|
|
637 r = r1 / (r1 + randg (b, 1))\n\ |
|
5730
|
638 @end example\n\ |
|
6547
|
639 @item @code{Erlang (a, n)}\n\ |
|
5730
|
640 @example\n\ |
|
6547
|
641 r = a * randg (n)\n\ |
|
5730
|
642 @end example\n\ |
|
6547
|
643 @item @code{chisq (df)} for @code{df > 0}\n\ |
|
5730
|
644 @example\n\ |
|
6547
|
645 r = 2 * randg (df / 2)\n\ |
|
5730
|
646 @end example\n\ |
|
|
647 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ |
|
|
648 @example\n\ |
|
6547
|
649 r = randn () / sqrt (2 * randg (df / 2) / df)\n\ |
|
5730
|
650 @end example\n\ |
|
6547
|
651 @item @code{F (n1, n2)} for @code{0 < n1}, @code{0 < n2}\n\ |
|
5730
|
652 @example\n\ |
|
7096
|
653 @group\n\ |
|
|
654 ## r1 equals 1 if n1 is infinite\n\ |
|
|
655 r1 = 2 * randg (n1 / 2) / n1\n\ |
|
|
656 ## r2 equals 1 if n2 is infinite\n\ |
|
|
657 r2 = 2 * randg (n2 / 2) / n2\n\ |
|
5730
|
658 r = r1 / r2\n\n\ |
|
7096
|
659 @end group\n\ |
|
5730
|
660 @end example\n\ |
|
|
661 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ |
|
|
662 @example\n\ |
|
6547
|
663 r = randp ((1 - p) / p * randg (n))\n\ |
|
5730
|
664 @end example\n\ |
|
6547
|
665 @item non-central @code{chisq (df, L)}, for @code{df >= 0} and @code{L > 0}\n\ |
|
5730
|
666 (use chisq if @code{L = 0})\n\ |
|
|
667 @example\n\ |
|
6547
|
668 r = randp (L / 2)\n\ |
|
|
669 r(r > 0) = 2 * randg (r(r > 0))\n\ |
|
|
670 r(df > 0) += 2 * randg (df(df > 0)/2)\n\ |
|
5730
|
671 @end example\n\ |
|
6547
|
672 @item @code{Dirichlet (a1, ..., ak)}\n\ |
|
5730
|
673 @example\n\ |
|
6547
|
674 r = (randg (a1), ..., randg (ak))\n\ |
|
|
675 r = r / sum (r)\n\ |
|
5730
|
676 @end example\n\ |
|
|
677 @end table\n\ |
|
6547
|
678 @seealso{rand, randn, rande, randp}\n\ |
|
5730
|
679 @end deftypefn") |
|
|
680 { |
|
|
681 octave_value retval; |
|
|
682 |
|
|
683 int nargin = args.length (); |
|
|
684 |
|
|
685 if (nargin < 1) |
|
|
686 error ("randg: insufficient arguments"); |
|
|
687 else |
|
|
688 { |
|
|
689 unwind_protect::begin_frame ("randg"); |
|
|
690 |
|
|
691 // This relies on the fact that elements are popped from the unwind |
|
|
692 // stack in the reverse of the order they are pushed |
|
|
693 // (i.e. current_distribution will be reset before calling |
|
|
694 // reset_rand_generator()). |
|
|
695 |
|
|
696 unwind_protect::add (reset_rand_generator, 0); |
|
|
697 unwind_protect_str (current_distribution); |
|
|
698 |
|
|
699 current_distribution = "gamma"; |
|
|
700 |
|
|
701 octave_rand::distribution (current_distribution); |
|
|
702 |
|
|
703 retval = do_rand (args, nargin, "randg", true); |
|
|
704 |
|
|
705 unwind_protect::run_frame ("randg"); |
|
|
706 } |
|
|
707 |
|
|
708 return retval; |
|
|
709 } |
|
|
710 |
|
|
711 /* |
|
|
712 %!test |
|
6437
|
713 %! randg("state",12) |
|
|
714 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
715 |
|
|
716 |
|
|
717 %!test |
|
|
718 %! % Test fixed state |
|
|
719 %! randg("state",1); |
|
6443
|
720 %! assert (randg(0.1,1,6), [0.0103951513331241 8.335671459898252e-05 0.00138691397249762 0.000587308416993855 0.495590518784736 2.3921917414795e-12],1e-6); |
|
6437
|
721 %!test |
|
|
722 %! % Test fixed state |
|
|
723 %! randg("state",1); |
|
6443
|
724 %! assert (randg(0.95,1,6), [3.099382433255327 0.3974529788871218 0.644367450750855 1.143261091802246 1.964111762696822 0.04011915547957939],1e-6); |
|
6437
|
725 %!test |
|
|
726 %! % Test fixed state |
|
|
727 %! randg("state",1); |
|
6443
|
728 %! assert (randg(1,1,6), [0.2273389379645993 1.288822625058359 0.2406335209340746 1.218869553370733 1.024649860162554 0.09631230343599533],1e-6); |
|
6437
|
729 %!test |
|
|
730 %! % Test fixed state |
|
|
731 %! randg("state",1); |
|
6443
|
732 %! assert (randg(10,1,6), [3.520369644331133 15.15369864472106 8.332112081991205 8.406211067432674 11.81193475187611 10.88792728177059],1e-5); |
|
6437
|
733 %!test |
|
|
734 %! % Test fixed state |
|
|
735 %! randg("state",1); |
|
6443
|
736 %! assert (randg(100,1,6), [75.34570255262264 115.4911985594699 95.23493031356388 95.48926019250911 106.2397448229803 103.4813150404118],1e-4); |
|
6437
|
737 %!test |
|
|
738 %! % Test fixed seed |
|
|
739 %! randg("seed",1); |
|
6443
|
740 %! assert (randg(0.1,1,6), [0.07144210487604141 0.460641473531723 0.4749028384685516 0.06823389977216721 0.000293838675133884 1.802567535340305e-12],1e-6); |
|
6437
|
741 %!test |
|
|
742 %! % Test fixed seed |
|
|
743 %! randg("seed",1); |
|
6443
|
744 %! assert (randg(0.95,1,6), [1.664905071258545 1.879976987838745 1.905677795410156 0.9948706030845642 0.5606933236122131 0.0766092911362648],1e-6); |
|
6437
|
745 %!test |
|
|
746 %! % Test fixed seed |
|
|
747 %! randg("seed",1); |
|
6443
|
748 %! assert (randg(1,1,6), [0.03512085229158401 0.6488978862762451 0.8114678859710693 0.1666885763406754 1.60791552066803 1.90356981754303],1e-6); |
|
6437
|
749 %!test |
|
|
750 %! % Test fixed seed |
|
|
751 %! randg("seed",1); |
|
6443
|
752 %! assert (randg(10,1,6), [6.566435813903809 10.11648464202881 10.73162078857422 7.747178077697754 6.278522491455078 6.240195751190186],1e-5); |
|
6437
|
753 %!test |
|
|
754 %! % Test fixed seed |
|
|
755 %! randg("seed",1); |
|
6443
|
756 %! assert (randg(100,1,6), [89.40208435058594 101.4734725952148 103.4020004272461 93.62763214111328 88.33104705810547 88.1871337890625],1e-4); |
|
6437
|
757 %!test |
|
|
758 %! if (__random_statistical_tests__) |
|
|
759 %! % statistical tests may fail occasionally. |
|
|
760 %! randg("state",12) |
|
|
761 %! a=0.1; x = randg(a,100000,1); |
|
|
762 %! assert(mean(x), a, 0.01); |
|
|
763 %! assert(var(x), a, 0.01); |
|
|
764 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
765 %! assert(kurtosis(x),6/a, 50.); |
|
|
766 %! endif |
|
|
767 %!test |
|
|
768 %! if (__random_statistical_tests__) |
|
|
769 %! % statistical tests may fail occasionally. |
|
|
770 %! randg("state",12) |
|
|
771 %! a=0.95; x = randg(a,100000,1); |
|
|
772 %! assert(mean(x), a, 0.01); |
|
|
773 %! assert(var(x), a, 0.04); |
|
|
774 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
775 %! assert(kurtosis(x),6/a, 2.); |
|
|
776 %! endif |
|
|
777 %!test |
|
|
778 %! if (__random_statistical_tests__) |
|
|
779 %! % statistical tests may fail occasionally. |
|
|
780 %! randg("state",12) |
|
|
781 %! a=1; x = randg(a,100000,1); |
|
|
782 %! assert(mean(x),a, 0.01); |
|
|
783 %! assert(var(x),a, 0.04); |
|
|
784 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
785 %! assert(kurtosis(x),6/a, 2.); |
|
|
786 %! endif |
|
|
787 %!test |
|
|
788 %! if (__random_statistical_tests__) |
|
|
789 %! % statistical tests may fail occasionally. |
|
|
790 %! randg("state",12) |
|
|
791 %! a=10; x = randg(a,100000,1); |
|
|
792 %! assert(mean(x), a, 0.1); |
|
|
793 %! assert(var(x), a, 0.5); |
|
|
794 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
795 %! assert(kurtosis(x),6/a, 0.5); |
|
|
796 %! endif |
|
|
797 %!test |
|
|
798 %! if (__random_statistical_tests__) |
|
|
799 %! % statistical tests may fail occasionally. |
|
|
800 %! randg("state",12) |
|
|
801 %! a=100; x = randg(a,100000,1); |
|
|
802 %! assert(mean(x), a, 0.2); |
|
|
803 %! assert(var(x), a, 2.); |
|
|
804 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
805 %! assert(kurtosis(x),6/a, 0.2); |
|
|
806 %! endif |
|
|
807 %!test |
|
|
808 %! randg("seed",12) |
|
5730
|
809 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
810 %!test |
|
6437
|
811 %! if (__random_statistical_tests__) |
|
|
812 %! % statistical tests may fail occasionally. |
|
|
813 %! randg("seed",12) |
|
|
814 %! a=0.1; x = randg(a,100000,1); |
|
|
815 %! assert(mean(x), a, 0.01); |
|
|
816 %! assert(var(x), a, 0.01); |
|
|
817 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
818 %! assert(kurtosis(x),6/a, 50.); |
|
|
819 %! endif |
|
5730
|
820 %!test |
|
6437
|
821 %! if (__random_statistical_tests__) |
|
|
822 %! % statistical tests may fail occasionally. |
|
|
823 %! randg("seed",12) |
|
|
824 %! a=0.95; x = randg(a,100000,1); |
|
|
825 %! assert(mean(x), a, 0.01); |
|
|
826 %! assert(var(x), a, 0.04); |
|
|
827 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
828 %! assert(kurtosis(x),6/a, 2.); |
|
|
829 %! endif |
|
5730
|
830 %!test |
|
6437
|
831 %! if (__random_statistical_tests__) |
|
|
832 %! % statistical tests may fail occasionally. |
|
|
833 %! randg("seed",12) |
|
|
834 %! a=1; x = randg(a,100000,1); |
|
|
835 %! assert(mean(x),a, 0.01); |
|
|
836 %! assert(var(x),a, 0.04); |
|
|
837 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
838 %! assert(kurtosis(x),6/a, 2.); |
|
|
839 %! endif |
|
5730
|
840 %!test |
|
6437
|
841 %! if (__random_statistical_tests__) |
|
|
842 %! % statistical tests may fail occasionally. |
|
|
843 %! randg("seed",12) |
|
|
844 %! a=10; x = randg(a,100000,1); |
|
|
845 %! assert(mean(x), a, 0.1); |
|
|
846 %! assert(var(x), a, 0.5); |
|
|
847 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
848 %! assert(kurtosis(x),6/a, 0.5); |
|
|
849 %! endif |
|
5730
|
850 %!test |
|
6437
|
851 %! if (__random_statistical_tests__) |
|
|
852 %! % statistical tests may fail occasionally. |
|
|
853 %! randg("seed",12) |
|
|
854 %! a=100; x = randg(a,100000,1); |
|
|
855 %! assert(mean(x), a, 0.2); |
|
|
856 %! assert(var(x), a, 2.); |
|
|
857 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
858 %! assert(kurtosis(x),6/a, 0.2); |
|
|
859 %! endif |
|
5730
|
860 */ |
|
|
861 |
|
|
862 |
|
|
863 DEFUN_DLD (randp, args, , |
|
|
864 "-*- texinfo -*-\n\ |
|
|
865 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ |
|
|
866 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ |
|
|
867 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ |
|
|
868 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ |
|
|
869 Return a matrix with Poisson distributed random elements. The arguments\n\ |
|
|
870 are handled the same as the arguments for @code{rand}, except for the\n\ |
|
|
871 argument @var{l}.\n\ |
|
|
872 \n\ |
|
|
873 Five different algorithms are used depending on the range of @var{l}\n\ |
|
|
874 and whether or not @var{l} is a scalar or a matrix.\n\ |
|
|
875 \n\ |
|
|
876 @table @asis\n\ |
|
|
877 @item For scalar @var{l} <= 12, use direct method.\n\ |
|
|
878 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
879 @item For scalar @var{l} > 12, use rejection method.[1]\n\ |
|
|
880 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
881 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ |
|
|
882 Stadlober E., et al., WinRand source code, available via FTP.\n\ |
|
|
883 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ |
|
|
884 Stadlober E., et al., WinRand source code, available via FTP, or\n\ |
|
|
885 H. Zechner, 'Efficient sampling from continuous and discrete\n\ |
|
|
886 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ |
|
|
887 University Graz, Austria, 1994.\n\ |
|
|
888 @item For @var{l} > 1e8, use normal approximation.\n\ |
|
|
889 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ |
|
|
890 D 50 p1284, 1994\n\ |
|
|
891 @end table\n\ |
|
6547
|
892 @seealso{rand, randn, rande, randg}\n\ |
|
5730
|
893 @end deftypefn") |
|
|
894 { |
|
|
895 octave_value retval; |
|
|
896 |
|
|
897 int nargin = args.length (); |
|
|
898 |
|
|
899 if (nargin < 1) |
|
|
900 error ("randp: insufficient arguments"); |
|
|
901 else |
|
|
902 { |
|
|
903 unwind_protect::begin_frame ("randp"); |
|
|
904 |
|
|
905 // This relies on the fact that elements are popped from the unwind |
|
|
906 // stack in the reverse of the order they are pushed |
|
|
907 // (i.e. current_distribution will be reset before calling |
|
|
908 // reset_rand_generator()). |
|
|
909 |
|
|
910 unwind_protect::add (reset_rand_generator, 0); |
|
|
911 unwind_protect_str (current_distribution); |
|
|
912 |
|
|
913 current_distribution = "poisson"; |
|
|
914 |
|
|
915 octave_rand::distribution (current_distribution); |
|
|
916 |
|
|
917 retval = do_rand (args, nargin, "randp", true); |
|
|
918 |
|
|
919 unwind_protect::run_frame ("randp"); |
|
|
920 } |
|
|
921 |
|
|
922 return retval; |
|
|
923 } |
|
|
924 |
|
|
925 /* |
|
|
926 %!test |
|
6437
|
927 %! randp("state",12) |
|
|
928 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
929 %!test |
|
|
930 %! % Test fixed state |
|
|
931 %! randp("state",1); |
|
|
932 %! assert(randp(5,1,6),[5 5 3 7 7 3]) |
|
|
933 %!test |
|
|
934 %! % Test fixed state |
|
|
935 %! randp("state",1); |
|
|
936 %! assert(randp(15,1,6),[13 15 8 18 18 15]) |
|
|
937 %!test |
|
|
938 %! % Test fixed state |
|
|
939 %! randp("state",1); |
|
6541
|
940 %! assert(randp(1e9,1,6),[999915677 999976657 1000047684 1000019035 999985749 999977692],1e-6) |
|
6437
|
941 %!test |
|
|
942 %! % Test fixed state |
|
|
943 %! randp("seed",1); |
|
6449
|
944 %! %%assert(randp(5,1,6),[8 2 3 6 6 8]) |
|
|
945 %! assert(randp(5,1,5),[8 2 3 6 6]) |
|
6437
|
946 %!test |
|
|
947 %! % Test fixed state |
|
|
948 %! randp("seed",1); |
|
|
949 %! assert(randp(15,1,6),[15 16 12 10 10 12]) |
|
|
950 %!test |
|
|
951 %! % Test fixed state |
|
|
952 %! randp("seed",1); |
|
|
953 %! assert(randp(1e9,1,6),[1000006208 1000012224 999981120 999963520 999963072 999981440]) |
|
|
954 %!test |
|
|
955 %! if (__random_statistical_tests__) |
|
|
956 %! % statistical tests may fail occasionally. |
|
|
957 %! randp("state",12) |
|
|
958 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
959 %! x = randp(a(1),100000,1); |
|
|
960 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
961 %! assert(mean(x),a(1),a(2)); |
|
|
962 %! assert(var(x),a(1),0.02*a(1)); |
|
|
963 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
964 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
965 %! endfor |
|
|
966 %! endif |
|
|
967 %!test |
|
|
968 %! if (__random_statistical_tests__) |
|
|
969 %! % statistical tests may fail occasionally. |
|
|
970 %! randp("state",12) |
|
|
971 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
972 %! x = randp(a(1)*ones(100000,1),100000,1); |
|
|
973 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
974 %! assert(mean(x),a(1),a(2)); |
|
|
975 %! assert(var(x),a(1),0.02*a(1)); |
|
|
976 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
977 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
978 %! endfor |
|
|
979 %! endif |
|
|
980 %!test |
|
|
981 %! randp("seed",12) |
|
5730
|
982 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
983 %!test |
|
6449
|
984 %! if (__random_statistical_tests__) |
|
|
985 %! % statistical tests may fail occasionally. |
|
|
986 %! randp("seed",12) |
|
|
987 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
988 %! x = randp(a(1),100000,1); |
|
|
989 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
990 %! assert(mean(x),a(1),a(2)); |
|
|
991 %! assert(var(x),a(1),0.02*a(1)); |
|
|
992 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
993 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
994 %! endfor |
|
|
995 %! endif |
|
5730
|
996 %!test |
|
6449
|
997 %! if (__random_statistical_tests__) |
|
|
998 %! % statistical tests may fail occasionally. |
|
|
999 %! randp("seed",12) |
|
|
1000 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
1001 %! x = randp(a(1)*ones(100000,1),100000,1); |
|
|
1002 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
1003 %! assert(mean(x),a(1),a(2)); |
|
|
1004 %! assert(var(x),a(1),0.02*a(1)); |
|
|
1005 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
1006 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
1007 %! endfor |
|
|
1008 %! endif |
|
5730
|
1009 */ |
|
|
1010 |
|
|
1011 /* |
|
2928
|
1012 ;;; Local Variables: *** |
|
|
1013 ;;; mode: C++ *** |
|
|
1014 ;;; End: *** |
|
|
1015 */ |
