|
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 |
|
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\ |
|
|
344 @code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\ |
|
|
345 (See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\ |
|
|
346 equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\ |
|
|
347 Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\ |
|
|
348 @url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\ |
|
6547
|
349 Do @strong{not} use for cryptography without securely hashing\n\ |
|
|
350 several returned values together, otherwise the generator state\n\ |
|
|
351 can be learned after reading 624 consecutive values.\n\ |
|
5730
|
352 \n\ |
|
|
353 @code{rand} includes a second random number generator, that was the\n\ |
|
5798
|
354 previous generator used in Octave. The new generator is used by default\n\ |
|
5730
|
355 as it is significantly faster than the old generator, and produces\n\ |
|
5798
|
356 random numbers with a significantly longer cycle time. However, in\n\ |
|
|
357 some circumstances it might be desirable to obtain the same random\n\ |
|
5730
|
358 sequences as used by the old generators. To do this the keyword\n\ |
|
|
359 \"seed\" is used to specify that the old generators should be use,\n\ |
|
|
360 as in\n\ |
|
2928
|
361 \n\ |
|
3369
|
362 @example\n\ |
|
5730
|
363 rand (\"seed\", val)\n\ |
|
3369
|
364 @end example\n\ |
|
|
365 \n\ |
|
5730
|
366 which sets the seed of the generator to @var{val}. The seed of the\n\ |
|
|
367 generator can be queried with\n\ |
|
|
368 \n\ |
|
|
369 @example\n\ |
|
|
370 s = rand (\"seed\")\n\ |
|
|
371 @end example\n\ |
|
|
372 \n\ |
|
|
373 However, it should be noted that querying the seed will not cause\n\ |
|
|
374 @code{rand} to use the old generators, only setting the seed will.\n\ |
|
|
375 To cause @code{rand} to once again use the new generators, the\n\ |
|
|
376 keyword \"state\" should be used to reset the state of the @code{rand}.\n\ |
|
5798
|
377 @seealso{randn, rande, randg, randp}\n\ |
|
3369
|
378 @end deftypefn") |
|
2928
|
379 { |
|
4307
|
380 octave_value retval; |
|
2928
|
381 |
|
|
382 int nargin = args.length (); |
|
|
383 |
|
4543
|
384 retval = do_rand (args, nargin, "rand"); |
|
2928
|
385 |
|
|
386 return retval; |
|
|
387 } |
|
|
388 |
|
5730
|
389 /* |
|
|
390 %!test # 'state' can be a scalar |
|
|
391 %! rand('state',12); x = rand(1,4); |
|
|
392 %! rand('state',12); y = rand(1,4); |
|
|
393 %! assert(x,y); |
|
|
394 %!test # 'state' can be a vector |
|
|
395 %! rand('state',[12,13]); x=rand(1,4); |
|
|
396 %! rand('state',[12;13]); y=rand(1,4); |
|
|
397 %! assert(x,y); |
|
|
398 %!test # querying 'state' doesn't disturb sequence |
|
|
399 %! rand('state',12); rand(1,2); x=rand(1,2); |
|
|
400 %! rand('state',12); rand(1,2); |
|
|
401 %! s=rand('state'); y=rand(1,2); |
|
|
402 %! assert(x,y); |
|
|
403 %! rand('state',s); z=rand(1,2); |
|
|
404 %! assert(x,z); |
|
|
405 %!test # 'seed' must be a scalar |
|
|
406 %! rand('seed',12); x = rand(1,4); |
|
|
407 %! rand('seed',12); y = rand(1,4); |
|
|
408 %! assert(x,y); |
|
|
409 %!error(rand('seed',[12,13])) |
|
|
410 %!test # querying 'seed' returns a value which can be used later |
|
|
411 %! s=rand('seed'); x=rand(1,2); |
|
|
412 %! rand('seed',s); y=rand(1,2); |
|
|
413 %! assert(x,y); |
|
|
414 %!test # querying 'seed' doesn't disturb sequence |
|
|
415 %! rand('seed',12); rand(1,2); x=rand(1,2); |
|
|
416 %! rand('seed',12); rand(1,2); |
|
|
417 %! s=rand('seed'); y=rand(1,2); |
|
|
418 %! assert(x,y); |
|
|
419 %! rand('seed',s); z=rand(1,2); |
|
|
420 %! assert(x,z); |
|
|
421 */ |
|
|
422 |
|
|
423 /* |
|
|
424 %!test |
|
6437
|
425 %! % Test fixed state |
|
|
426 %! rand("state",1); |
|
6443
|
427 %! assert (rand(1,6), [0.1343642441124013 0.8474337369372327 0.763774618976614 0.2550690257394218 0.495435087091941 0.4494910647887382],1e-6); |
|
6437
|
428 %!test |
|
6443
|
429 %! % Test fixed seed |
|
6437
|
430 %! rand("seed",1); |
|
6443
|
431 %! assert (rand(1,6), [0.8668024251237512 0.9126510815694928 0.09366085007786751 0.1664607301354408 0.7408077004365623 0.7615650338120759],1e-6); |
|
5730
|
432 %!test |
|
6437
|
433 %! if (__random_statistical_tests__) |
|
|
434 %! % statistical tests may fail occasionally. |
|
|
435 %! rand("state",12); |
|
|
436 %! x = rand(100000,1); |
|
|
437 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
438 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
439 %! assert(mean(x),0.5,0.0024); |
|
|
440 %! assert(var(x),1/48,0.0632); |
|
|
441 %! assert(skewness(x),0,0.012); |
|
|
442 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
443 %! endif |
|
|
444 %!test |
|
|
445 %! if (__random_statistical_tests__) |
|
|
446 %! % statistical tests may fail occasionally. |
|
|
447 %! rand("seed",12); |
|
|
448 %! x = rand(100000,1); |
|
|
449 %! assert(max(x)<1.); %*** Please report this!!! *** |
|
|
450 %! assert(min(x)>0.); %*** Please report this!!! *** |
|
|
451 %! assert(mean(x),0.5,0.0024); |
|
|
452 %! assert(var(x),1/48,0.0632); |
|
|
453 %! assert(skewness(x),0,0.012); |
|
|
454 %! assert(kurtosis(x),-6/5,0.0094); |
|
|
455 %! endif |
|
5730
|
456 */ |
|
|
457 |
|
|
458 |
|
4307
|
459 static std::string current_distribution = octave_rand::distribution (); |
|
|
460 |
|
2928
|
461 static void |
|
|
462 reset_rand_generator (void *) |
|
|
463 { |
|
4307
|
464 octave_rand::distribution (current_distribution); |
|
2928
|
465 } |
|
|
466 |
|
4665
|
467 DEFUN_DLD (randn, args, , |
|
3369
|
468 "-*- texinfo -*-\n\ |
|
|
469 @deftypefn {Loadable Function} {} randn (@var{x})\n\ |
|
|
470 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\ |
|
5730
|
471 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\ |
|
|
472 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\ |
|
|
473 Return a matrix with normally distributed random elements. The\n\ |
|
|
474 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
3369
|
475 \n\ |
|
5730
|
476 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
477 transform from a uniform to a normal distribution. (G. Marsaglia and\n\ |
|
|
478 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
479 J. Statistical Software, vol 5, 2000,\n\ |
|
|
480 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
2928
|
481 \n\ |
|
6547
|
482 @seealso{rand, rande, randg, randp}\n\ |
|
3369
|
483 @end deftypefn") |
|
2928
|
484 { |
|
4307
|
485 octave_value retval; |
|
2928
|
486 |
|
|
487 int nargin = args.length (); |
|
|
488 |
|
4543
|
489 unwind_protect::begin_frame ("randn"); |
|
2928
|
490 |
|
4543
|
491 // This relies on the fact that elements are popped from the unwind |
|
|
492 // stack in the reverse of the order they are pushed |
|
|
493 // (i.e. current_distribution will be reset before calling |
|
|
494 // reset_rand_generator()). |
|
2928
|
495 |
|
4543
|
496 unwind_protect::add (reset_rand_generator, 0); |
|
|
497 unwind_protect_str (current_distribution); |
|
2928
|
498 |
|
4543
|
499 current_distribution = "normal"; |
|
2928
|
500 |
|
4543
|
501 octave_rand::distribution (current_distribution); |
|
2928
|
502 |
|
4543
|
503 retval = do_rand (args, nargin, "randn"); |
|
2928
|
504 |
|
4543
|
505 unwind_protect::run_frame ("randn"); |
|
2928
|
506 |
|
|
507 return retval; |
|
|
508 } |
|
|
509 |
|
|
510 /* |
|
5730
|
511 %!test |
|
6437
|
512 %! % Test fixed state |
|
|
513 %! randn("state",1); |
|
6443
|
514 %! assert (randn(1,6), [-2.666521678978671 -0.7381719971724564 1.507903992673601 0.6019427189162239 -0.450661261143348 -0.7054431351574116],1e-6); |
|
6437
|
515 %!test |
|
6443
|
516 %! % Test fixed seed |
|
6437
|
517 %! randn("seed",1); |
|
6443
|
518 %! assert (randn(1,6), [-1.039402365684509 -1.25938892364502 0.1968704611063004 0.3874166905879974 -0.5976632833480835 -0.6615074276924133],1e-6); |
|
5730
|
519 %!test |
|
6437
|
520 %! if (__random_statistical_tests__) |
|
|
521 %! % statistical tests may fail occasionally. |
|
|
522 %! randn("state",12); |
|
|
523 %! x = randn(100000,1); |
|
|
524 %! assert(mean(x),0,0.01); |
|
|
525 %! assert(var(x),1,0.02); |
|
|
526 %! assert(skewness(x),0,0.02); |
|
|
527 %! assert(kurtosis(x),0,0.04); |
|
|
528 %! endif |
|
|
529 %!test |
|
|
530 %! if (__random_statistical_tests__) |
|
|
531 %! % statistical tests may fail occasionally. |
|
|
532 %! randn("seed",12); |
|
|
533 %! x = randn(100000,1); |
|
|
534 %! assert(mean(x),0,0.01); |
|
|
535 %! assert(var(x),1,0.02); |
|
|
536 %! assert(skewness(x),0,0.02); |
|
|
537 %! assert(kurtosis(x),0,0.04); |
|
|
538 %! endif |
|
5730
|
539 */ |
|
|
540 |
|
|
541 DEFUN_DLD (rande, args, , |
|
|
542 "-*- texinfo -*-\n\ |
|
|
543 @deftypefn {Loadable Function} {} rande (@var{x})\n\ |
|
|
544 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\ |
|
|
545 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\ |
|
|
546 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\ |
|
|
547 Return a matrix with exponentially distributed random elements. The\n\ |
|
|
548 arguments are handled the same as the arguments for @code{rand}.\n\ |
|
|
549 \n\ |
|
|
550 By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\ |
|
|
551 transform from a uniform to a exponential distribution. (G. Marsaglia and\n\ |
|
|
552 W.W. Tsang, 'Ziggurat method for generating random variables',\n\ |
|
|
553 J. Statistical Software, vol 5, 2000,\n\ |
|
|
554 @url{http://www.jstatsoft.org/v05/i08/})\n\ |
|
6547
|
555 @seealso{rand, randn, randg, randp}\n\ |
|
5730
|
556 @end deftypefn") |
|
|
557 { |
|
|
558 octave_value retval; |
|
|
559 |
|
|
560 int nargin = args.length (); |
|
|
561 |
|
|
562 unwind_protect::begin_frame ("rande"); |
|
|
563 |
|
|
564 // This relies on the fact that elements are popped from the unwind |
|
|
565 // stack in the reverse of the order they are pushed |
|
|
566 // (i.e. current_distribution will be reset before calling |
|
|
567 // reset_rand_generator()). |
|
|
568 |
|
|
569 unwind_protect::add (reset_rand_generator, 0); |
|
|
570 unwind_protect_str (current_distribution); |
|
|
571 |
|
|
572 current_distribution = "exponential"; |
|
|
573 |
|
|
574 octave_rand::distribution (current_distribution); |
|
|
575 |
|
|
576 retval = do_rand (args, nargin, "rande"); |
|
|
577 |
|
|
578 unwind_protect::run_frame ("rande"); |
|
|
579 |
|
|
580 return retval; |
|
|
581 } |
|
|
582 |
|
|
583 /* |
|
|
584 %!test |
|
6437
|
585 %! % Test fixed state |
|
|
586 %! rande("state",1); |
|
6443
|
587 %! assert (rande(1,6), [3.602973885835625 0.1386190677555021 0.6743112889616958 0.4512830847258422 0.7255744741233175 0.3415969205292291],1e-6); |
|
6437
|
588 %!test |
|
6443
|
589 %! % Test fixed seed |
|
6437
|
590 %! rande("seed",1); |
|
6443
|
591 %! assert (rande(1,6), [0.06492075175653866 1.717980206012726 0.4816154008731246 0.5231300676241517 0.103910739364359 1.668931916356087],1e-6); |
|
5730
|
592 %!test |
|
6437
|
593 %! if (__random_statistical_tests__) |
|
|
594 %! % statistical tests may fail occasionally |
|
|
595 %! rande("state",1); |
|
|
596 %! x = rande(100000,1); |
|
|
597 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
598 %! assert(mean(x),1,0.01); |
|
|
599 %! assert(var(x),1,0.03); |
|
|
600 %! assert(skewness(x),2,0.06); |
|
|
601 %! assert(kurtosis(x),6,0.7); |
|
|
602 %! endif |
|
|
603 %!test |
|
|
604 %! if (__random_statistical_tests__) |
|
|
605 %! % statistical tests may fail occasionally |
|
|
606 %! rande("seed",1); |
|
|
607 %! x = rande(100000,1); |
|
|
608 %! assert(min(x)>0); % *** Please report this!!! *** |
|
|
609 %! assert(mean(x),1,0.01); |
|
|
610 %! assert(var(x),1,0.03); |
|
|
611 %! assert(skewness(x),2,0.06); |
|
|
612 %! assert(kurtosis(x),6,0.7); |
|
|
613 %! endif |
|
5730
|
614 */ |
|
|
615 |
|
|
616 DEFUN_DLD (randg, args, , |
|
|
617 "-*- texinfo -*-\n\ |
|
|
618 @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\ |
|
|
619 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\ |
|
|
620 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\ |
|
|
621 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\ |
|
|
622 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\ |
|
|
623 The arguments are handled the same as the arguments for @code{rand},\n\ |
|
|
624 except for the argument @var{a}.\n\ |
|
|
625 \n\ |
|
|
626 This can be used to generate many distributions:\n\ |
|
|
627 \n\ |
|
|
628 @table @asis\n\ |
|
6547
|
629 @item @code{gamma (a, b)} for @code{a > -1}, @code{b > 0}\n\ |
|
5730
|
630 @example\n\ |
|
6547
|
631 r = b * randg (a)\n\ |
|
5730
|
632 @end example\n\ |
|
6547
|
633 @item @code{beta (a, b)} for @code{a > -1}, @code{b > -1}\n\ |
|
5730
|
634 @example\n\ |
|
6547
|
635 r1 = randg (a, 1)\n\ |
|
|
636 r = r1 / (r1 + randg (b, 1))\n\ |
|
5730
|
637 @end example\n\ |
|
6547
|
638 @item @code{Erlang (a, n)}\n\ |
|
5730
|
639 @example\n\ |
|
6547
|
640 r = a * randg (n)\n\ |
|
5730
|
641 @end example\n\ |
|
6547
|
642 @item @code{chisq (df)} for @code{df > 0}\n\ |
|
5730
|
643 @example\n\ |
|
6547
|
644 r = 2 * randg (df / 2)\n\ |
|
5730
|
645 @end example\n\ |
|
|
646 @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\ |
|
|
647 @example\n\ |
|
6547
|
648 r = randn () / sqrt (2 * randg (df / 2) / df)\n\ |
|
5730
|
649 @end example\n\ |
|
6547
|
650 @item @code{F (n1, n2)} for @code{0 < n1}, @code{0 < n2}\n\ |
|
5730
|
651 @example\n\ |
|
6547
|
652 r1 = 2 * randg (n1 / 2) / n1 or 1 if n1 is infinite\n\ |
|
|
653 r2 = 2 * randg (n2 / 2) / n2 or 1 if n2 is infinite\n\ |
|
5730
|
654 r = r1 / r2\n\n\ |
|
|
655 @end example\n\ |
|
|
656 @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\ |
|
|
657 @example\n\ |
|
6547
|
658 r = randp ((1 - p) / p * randg (n))\n\ |
|
5730
|
659 @end example\n\ |
|
6547
|
660 @item non-central @code{chisq (df, L)}, for @code{df >= 0} and @code{L > 0}\n\ |
|
5730
|
661 (use chisq if @code{L = 0})\n\ |
|
|
662 @example\n\ |
|
6547
|
663 r = randp (L / 2)\n\ |
|
|
664 r(r > 0) = 2 * randg (r(r > 0))\n\ |
|
|
665 r(df > 0) += 2 * randg (df(df > 0)/2)\n\ |
|
5730
|
666 @end example\n\ |
|
6547
|
667 @item @code{Dirichlet (a1, ..., ak)}\n\ |
|
5730
|
668 @example\n\ |
|
6547
|
669 r = (randg (a1), ..., randg (ak))\n\ |
|
|
670 r = r / sum (r)\n\ |
|
5730
|
671 @end example\n\ |
|
|
672 @end table\n\ |
|
6547
|
673 @seealso{rand, randn, rande, randp}\n\ |
|
5730
|
674 @end deftypefn") |
|
|
675 { |
|
|
676 octave_value retval; |
|
|
677 |
|
|
678 int nargin = args.length (); |
|
|
679 |
|
|
680 if (nargin < 1) |
|
|
681 error ("randg: insufficient arguments"); |
|
|
682 else |
|
|
683 { |
|
|
684 unwind_protect::begin_frame ("randg"); |
|
|
685 |
|
|
686 // This relies on the fact that elements are popped from the unwind |
|
|
687 // stack in the reverse of the order they are pushed |
|
|
688 // (i.e. current_distribution will be reset before calling |
|
|
689 // reset_rand_generator()). |
|
|
690 |
|
|
691 unwind_protect::add (reset_rand_generator, 0); |
|
|
692 unwind_protect_str (current_distribution); |
|
|
693 |
|
|
694 current_distribution = "gamma"; |
|
|
695 |
|
|
696 octave_rand::distribution (current_distribution); |
|
|
697 |
|
|
698 retval = do_rand (args, nargin, "randg", true); |
|
|
699 |
|
|
700 unwind_protect::run_frame ("randg"); |
|
|
701 } |
|
|
702 |
|
|
703 return retval; |
|
|
704 } |
|
|
705 |
|
|
706 /* |
|
|
707 %!test |
|
6437
|
708 %! randg("state",12) |
|
|
709 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
710 |
|
|
711 |
|
|
712 %!test |
|
|
713 %! % Test fixed state |
|
|
714 %! randg("state",1); |
|
6443
|
715 %! assert (randg(0.1,1,6), [0.0103951513331241 8.335671459898252e-05 0.00138691397249762 0.000587308416993855 0.495590518784736 2.3921917414795e-12],1e-6); |
|
6437
|
716 %!test |
|
|
717 %! % Test fixed state |
|
|
718 %! randg("state",1); |
|
6443
|
719 %! assert (randg(0.95,1,6), [3.099382433255327 0.3974529788871218 0.644367450750855 1.143261091802246 1.964111762696822 0.04011915547957939],1e-6); |
|
6437
|
720 %!test |
|
|
721 %! % Test fixed state |
|
|
722 %! randg("state",1); |
|
6443
|
723 %! assert (randg(1,1,6), [0.2273389379645993 1.288822625058359 0.2406335209340746 1.218869553370733 1.024649860162554 0.09631230343599533],1e-6); |
|
6437
|
724 %!test |
|
|
725 %! % Test fixed state |
|
|
726 %! randg("state",1); |
|
6443
|
727 %! assert (randg(10,1,6), [3.520369644331133 15.15369864472106 8.332112081991205 8.406211067432674 11.81193475187611 10.88792728177059],1e-5); |
|
6437
|
728 %!test |
|
|
729 %! % Test fixed state |
|
|
730 %! randg("state",1); |
|
6443
|
731 %! assert (randg(100,1,6), [75.34570255262264 115.4911985594699 95.23493031356388 95.48926019250911 106.2397448229803 103.4813150404118],1e-4); |
|
6437
|
732 %!test |
|
|
733 %! % Test fixed seed |
|
|
734 %! randg("seed",1); |
|
6443
|
735 %! assert (randg(0.1,1,6), [0.07144210487604141 0.460641473531723 0.4749028384685516 0.06823389977216721 0.000293838675133884 1.802567535340305e-12],1e-6); |
|
6437
|
736 %!test |
|
|
737 %! % Test fixed seed |
|
|
738 %! randg("seed",1); |
|
6443
|
739 %! assert (randg(0.95,1,6), [1.664905071258545 1.879976987838745 1.905677795410156 0.9948706030845642 0.5606933236122131 0.0766092911362648],1e-6); |
|
6437
|
740 %!test |
|
|
741 %! % Test fixed seed |
|
|
742 %! randg("seed",1); |
|
6443
|
743 %! assert (randg(1,1,6), [0.03512085229158401 0.6488978862762451 0.8114678859710693 0.1666885763406754 1.60791552066803 1.90356981754303],1e-6); |
|
6437
|
744 %!test |
|
|
745 %! % Test fixed seed |
|
|
746 %! randg("seed",1); |
|
6443
|
747 %! assert (randg(10,1,6), [6.566435813903809 10.11648464202881 10.73162078857422 7.747178077697754 6.278522491455078 6.240195751190186],1e-5); |
|
6437
|
748 %!test |
|
|
749 %! % Test fixed seed |
|
|
750 %! randg("seed",1); |
|
6443
|
751 %! assert (randg(100,1,6), [89.40208435058594 101.4734725952148 103.4020004272461 93.62763214111328 88.33104705810547 88.1871337890625],1e-4); |
|
6437
|
752 %!test |
|
|
753 %! if (__random_statistical_tests__) |
|
|
754 %! % statistical tests may fail occasionally. |
|
|
755 %! randg("state",12) |
|
|
756 %! a=0.1; x = randg(a,100000,1); |
|
|
757 %! assert(mean(x), a, 0.01); |
|
|
758 %! assert(var(x), a, 0.01); |
|
|
759 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
760 %! assert(kurtosis(x),6/a, 50.); |
|
|
761 %! endif |
|
|
762 %!test |
|
|
763 %! if (__random_statistical_tests__) |
|
|
764 %! % statistical tests may fail occasionally. |
|
|
765 %! randg("state",12) |
|
|
766 %! a=0.95; x = randg(a,100000,1); |
|
|
767 %! assert(mean(x), a, 0.01); |
|
|
768 %! assert(var(x), a, 0.04); |
|
|
769 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
770 %! assert(kurtosis(x),6/a, 2.); |
|
|
771 %! endif |
|
|
772 %!test |
|
|
773 %! if (__random_statistical_tests__) |
|
|
774 %! % statistical tests may fail occasionally. |
|
|
775 %! randg("state",12) |
|
|
776 %! a=1; x = randg(a,100000,1); |
|
|
777 %! assert(mean(x),a, 0.01); |
|
|
778 %! assert(var(x),a, 0.04); |
|
|
779 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
780 %! assert(kurtosis(x),6/a, 2.); |
|
|
781 %! endif |
|
|
782 %!test |
|
|
783 %! if (__random_statistical_tests__) |
|
|
784 %! % statistical tests may fail occasionally. |
|
|
785 %! randg("state",12) |
|
|
786 %! a=10; x = randg(a,100000,1); |
|
|
787 %! assert(mean(x), a, 0.1); |
|
|
788 %! assert(var(x), a, 0.5); |
|
|
789 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
790 %! assert(kurtosis(x),6/a, 0.5); |
|
|
791 %! endif |
|
|
792 %!test |
|
|
793 %! if (__random_statistical_tests__) |
|
|
794 %! % statistical tests may fail occasionally. |
|
|
795 %! randg("state",12) |
|
|
796 %! a=100; x = randg(a,100000,1); |
|
|
797 %! assert(mean(x), a, 0.2); |
|
|
798 %! assert(var(x), a, 2.); |
|
|
799 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
800 %! assert(kurtosis(x),6/a, 0.2); |
|
|
801 %! endif |
|
|
802 %!test |
|
|
803 %! randg("seed",12) |
|
5730
|
804 %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report |
|
|
805 %!test |
|
6437
|
806 %! if (__random_statistical_tests__) |
|
|
807 %! % statistical tests may fail occasionally. |
|
|
808 %! randg("seed",12) |
|
|
809 %! a=0.1; x = randg(a,100000,1); |
|
|
810 %! assert(mean(x), a, 0.01); |
|
|
811 %! assert(var(x), a, 0.01); |
|
|
812 %! assert(skewness(x),2/sqrt(a), 1.); |
|
|
813 %! assert(kurtosis(x),6/a, 50.); |
|
|
814 %! endif |
|
5730
|
815 %!test |
|
6437
|
816 %! if (__random_statistical_tests__) |
|
|
817 %! % statistical tests may fail occasionally. |
|
|
818 %! randg("seed",12) |
|
|
819 %! a=0.95; x = randg(a,100000,1); |
|
|
820 %! assert(mean(x), a, 0.01); |
|
|
821 %! assert(var(x), a, 0.04); |
|
|
822 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
823 %! assert(kurtosis(x),6/a, 2.); |
|
|
824 %! endif |
|
5730
|
825 %!test |
|
6437
|
826 %! if (__random_statistical_tests__) |
|
|
827 %! % statistical tests may fail occasionally. |
|
|
828 %! randg("seed",12) |
|
|
829 %! a=1; x = randg(a,100000,1); |
|
|
830 %! assert(mean(x),a, 0.01); |
|
|
831 %! assert(var(x),a, 0.04); |
|
|
832 %! assert(skewness(x),2/sqrt(a), 0.2); |
|
|
833 %! assert(kurtosis(x),6/a, 2.); |
|
|
834 %! endif |
|
5730
|
835 %!test |
|
6437
|
836 %! if (__random_statistical_tests__) |
|
|
837 %! % statistical tests may fail occasionally. |
|
|
838 %! randg("seed",12) |
|
|
839 %! a=10; x = randg(a,100000,1); |
|
|
840 %! assert(mean(x), a, 0.1); |
|
|
841 %! assert(var(x), a, 0.5); |
|
|
842 %! assert(skewness(x),2/sqrt(a), 0.1); |
|
|
843 %! assert(kurtosis(x),6/a, 0.5); |
|
|
844 %! endif |
|
5730
|
845 %!test |
|
6437
|
846 %! if (__random_statistical_tests__) |
|
|
847 %! % statistical tests may fail occasionally. |
|
|
848 %! randg("seed",12) |
|
|
849 %! a=100; x = randg(a,100000,1); |
|
|
850 %! assert(mean(x), a, 0.2); |
|
|
851 %! assert(var(x), a, 2.); |
|
|
852 %! assert(skewness(x),2/sqrt(a), 0.05); |
|
|
853 %! assert(kurtosis(x),6/a, 0.2); |
|
|
854 %! endif |
|
5730
|
855 */ |
|
|
856 |
|
|
857 |
|
|
858 DEFUN_DLD (randp, args, , |
|
|
859 "-*- texinfo -*-\n\ |
|
|
860 @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\ |
|
|
861 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\ |
|
|
862 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\ |
|
|
863 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\ |
|
|
864 Return a matrix with Poisson distributed random elements. The arguments\n\ |
|
|
865 are handled the same as the arguments for @code{rand}, except for the\n\ |
|
|
866 argument @var{l}.\n\ |
|
|
867 \n\ |
|
|
868 Five different algorithms are used depending on the range of @var{l}\n\ |
|
|
869 and whether or not @var{l} is a scalar or a matrix.\n\ |
|
|
870 \n\ |
|
|
871 @table @asis\n\ |
|
|
872 @item For scalar @var{l} <= 12, use direct method.\n\ |
|
|
873 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
874 @item For scalar @var{l} > 12, use rejection method.[1]\n\ |
|
|
875 Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\ |
|
|
876 @item For matrix @var{l} <= 10, use inversion method.[2]\n\ |
|
|
877 Stadlober E., et al., WinRand source code, available via FTP.\n\ |
|
|
878 @item For matrix @var{l} > 10, use patchwork rejection method.\n\ |
|
|
879 Stadlober E., et al., WinRand source code, available via FTP, or\n\ |
|
|
880 H. Zechner, 'Efficient sampling from continuous and discrete\n\ |
|
|
881 unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\ |
|
|
882 University Graz, Austria, 1994.\n\ |
|
|
883 @item For @var{l} > 1e8, use normal approximation.\n\ |
|
|
884 L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\ |
|
|
885 D 50 p1284, 1994\n\ |
|
|
886 @end table\n\ |
|
6547
|
887 @seealso{rand, randn, rande, randg}\n\ |
|
5730
|
888 @end deftypefn") |
|
|
889 { |
|
|
890 octave_value retval; |
|
|
891 |
|
|
892 int nargin = args.length (); |
|
|
893 |
|
|
894 if (nargin < 1) |
|
|
895 error ("randp: insufficient arguments"); |
|
|
896 else |
|
|
897 { |
|
|
898 unwind_protect::begin_frame ("randp"); |
|
|
899 |
|
|
900 // This relies on the fact that elements are popped from the unwind |
|
|
901 // stack in the reverse of the order they are pushed |
|
|
902 // (i.e. current_distribution will be reset before calling |
|
|
903 // reset_rand_generator()). |
|
|
904 |
|
|
905 unwind_protect::add (reset_rand_generator, 0); |
|
|
906 unwind_protect_str (current_distribution); |
|
|
907 |
|
|
908 current_distribution = "poisson"; |
|
|
909 |
|
|
910 octave_rand::distribution (current_distribution); |
|
|
911 |
|
|
912 retval = do_rand (args, nargin, "randp", true); |
|
|
913 |
|
|
914 unwind_protect::run_frame ("randp"); |
|
|
915 } |
|
|
916 |
|
|
917 return retval; |
|
|
918 } |
|
|
919 |
|
|
920 /* |
|
|
921 %!test |
|
6437
|
922 %! randp("state",12) |
|
|
923 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
924 %!test |
|
|
925 %! % Test fixed state |
|
|
926 %! randp("state",1); |
|
|
927 %! assert(randp(5,1,6),[5 5 3 7 7 3]) |
|
|
928 %!test |
|
|
929 %! % Test fixed state |
|
|
930 %! randp("state",1); |
|
|
931 %! assert(randp(15,1,6),[13 15 8 18 18 15]) |
|
|
932 %!test |
|
|
933 %! % Test fixed state |
|
|
934 %! randp("state",1); |
|
6541
|
935 %! assert(randp(1e9,1,6),[999915677 999976657 1000047684 1000019035 999985749 999977692],1e-6) |
|
6437
|
936 %!test |
|
|
937 %! % Test fixed state |
|
|
938 %! randp("seed",1); |
|
6449
|
939 %! %%assert(randp(5,1,6),[8 2 3 6 6 8]) |
|
|
940 %! assert(randp(5,1,5),[8 2 3 6 6]) |
|
6437
|
941 %!test |
|
|
942 %! % Test fixed state |
|
|
943 %! randp("seed",1); |
|
|
944 %! assert(randp(15,1,6),[15 16 12 10 10 12]) |
|
|
945 %!test |
|
|
946 %! % Test fixed state |
|
|
947 %! randp("seed",1); |
|
|
948 %! assert(randp(1e9,1,6),[1000006208 1000012224 999981120 999963520 999963072 999981440]) |
|
|
949 %!test |
|
|
950 %! if (__random_statistical_tests__) |
|
|
951 %! % statistical tests may fail occasionally. |
|
|
952 %! randp("state",12) |
|
|
953 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
954 %! x = randp(a(1),100000,1); |
|
|
955 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
956 %! assert(mean(x),a(1),a(2)); |
|
|
957 %! assert(var(x),a(1),0.02*a(1)); |
|
|
958 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
959 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
960 %! endfor |
|
|
961 %! endif |
|
|
962 %!test |
|
|
963 %! if (__random_statistical_tests__) |
|
|
964 %! % statistical tests may fail occasionally. |
|
|
965 %! randp("state",12) |
|
|
966 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
967 %! x = randp(a(1)*ones(100000,1),100000,1); |
|
|
968 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
969 %! assert(mean(x),a(1),a(2)); |
|
|
970 %! assert(var(x),a(1),0.02*a(1)); |
|
|
971 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
972 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
973 %! endfor |
|
|
974 %! endif |
|
|
975 %!test |
|
|
976 %! randp("seed",12) |
|
5730
|
977 %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report |
|
|
978 %!test |
|
6449
|
979 %! if (__random_statistical_tests__) |
|
|
980 %! % statistical tests may fail occasionally. |
|
|
981 %! randp("seed",12) |
|
|
982 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
983 %! x = randp(a(1),100000,1); |
|
|
984 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
985 %! assert(mean(x),a(1),a(2)); |
|
|
986 %! assert(var(x),a(1),0.02*a(1)); |
|
|
987 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
988 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
989 %! endfor |
|
|
990 %! endif |
|
5730
|
991 %!test |
|
6449
|
992 %! if (__random_statistical_tests__) |
|
|
993 %! % statistical tests may fail occasionally. |
|
|
994 %! randp("seed",12) |
|
|
995 %! for a=[5, 15, 1e9; 0.03, 0.03, -5e-3; 0.03, 0.03, 0.03] |
|
|
996 %! x = randp(a(1)*ones(100000,1),100000,1); |
|
|
997 %! assert(min(x)>=0); % *** Please report this!!! *** |
|
|
998 %! assert(mean(x),a(1),a(2)); |
|
|
999 %! assert(var(x),a(1),0.02*a(1)); |
|
|
1000 %! assert(skewness(x),1/sqrt(a(1)),a(3)); |
|
|
1001 %! assert(kurtosis(x),1/a(1),3*a(3)); |
|
|
1002 %! endfor |
|
|
1003 %! endif |
|
5730
|
1004 */ |
|
|
1005 |
|
|
1006 /* |
|
2928
|
1007 ;;; Local Variables: *** |
|
|
1008 ;;; mode: C++ *** |
|
|
1009 ;;; End: *** |
|
|
1010 */ |
