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