~ubuntu-branches/ubuntu/utopic/libmath-random-perl/utopic : revision 5

1

package Math::Random;

2

3

use strict;

4

use Carp;

5

use vars qw($VERSION @ISA @EXPORT @EXPORT_OK %EXPORT_TAGS $AUTOLOAD);

6

7

require Exporter;

8

require DynaLoader;

9

require AutoLoader;

10

11

@ISA = qw(Exporter DynaLoader);

12

$VERSION = '0.71';

13

14

@EXPORT = qw(random_normal

15

random_permutation

16

random_permuted_index

17

random_uniform

18

random_uniform_integer

19

random_seed_from_phrase

20

random_get_seed

21

random_set_seed_from_phrase

22

random_set_seed

23

);

24

25

@EXPORT_OK = qw(random_beta

26

random_chi_square

27

random_exponential

28

random_f

29

random_gamma

30

random_multivariate_normal

31

random_multinomial

32

random_noncentral_chi_square

33

random_noncentral_f

34

random_normal

35

random_permutation

36

random_permuted_index

37

random_uniform

38

random_poisson

39

random_uniform_integer

40

random_negative_binomial

41

random_binomial

42

random_seed_from_phrase

43

random_get_seed

44

random_set_seed_from_phrase

45

random_set_seed

46

);

47

48

%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

49

50

sub AUTOLOAD {

51

# This AUTOLOAD is used to 'autoload' constants from the constant()

52

# XS function. If a constant is not found then control is passed

53

# to the AUTOLOAD in AutoLoader.

54

55

my $constname;

56

($constname = $AUTOLOAD) =~ s/.*:://;

57

croak "& not defined" if $constname eq 'constant';

58

my $val = constant($constname, @_ ? $_[0] : 0);

59

if ($! != 0) {

60

if ($! =~ /Invalid/) {

61

$AutoLoader::AUTOLOAD = $AUTOLOAD;

62

goto &AutoLoader::AUTOLOAD;

63

}

64

else {

65

croak "Your vendor has not defined Math::Random macro $constname";

66

}

67

}

68

*$AUTOLOAD = sub () { $val };

69

goto &$AUTOLOAD;

70

}

71

72

bootstrap Math::Random $VERSION;

73

74

75

### set seeds by default

76

salfph(scalar(localtime()));

77

78

#####################################################################

79

# RANDOM DEVIATE GENERATORS #

80

#####################################################################

81

82

sub random_beta { # Arguments: ($n,$aa,$bb)

83

croak "Usage: random_beta(\$n,\$aa,\$bb)" if scalar(@_) < 3;

84

my($n, $aa, $bb) = @_;

85

croak("($aa = \$aa < 1.0E-37) or ($bb = \$bb < 1.0E-37)\nin ".

86

"random_beta(\$n,\$aa,\$bb)")

87

if (($aa < 1.0E-37) or ($bb < 1.0E-37));

88

return genbet($aa,$bb) unless wantarray();

89

my $val;

90

my @ans = (0) x $n;

91

foreach $val (@ans) { $val = genbet($aa,$bb); }

92

return @ans;

93

}

94

95

sub random_chi_square { # Arguments: ($n,$df)

96

croak "Usage: random_chi_square(\$n,\$df)" if scalar(@_) < 2;

97

my($n, $df) = @_;

98

croak "$df = \$df <= 0\nin random_chi_square(\$n,\$df)" if ($df <= 0);

99

return genchi($df) unless wantarray();

100

my $val;

101

my @ans = (0) x $n;

102

foreach $val (@ans) { $val = genchi($df); }

103

return @ans;

104

}

105

106

sub random_exponential { # Arguments: ($n,$av), defaults (1,1)

107

return wantarray() ? (genexp(1)) : genexp(1)

108

if scalar(@_) == 0; # default behavior if no arguments

109

my($n, $av) = @_;

110

$av = 1 unless defined($av); # default $av is 1

111

croak "$av = \$av < 0\nin random_exponential(\$n,\$av)" if ($av < 0);

112

return genexp($av) unless wantarray();

113

my $val;

114

my @ans = (0) x $n;

115

foreach $val (@ans) { $val = genexp($av); }

116

return @ans;

117

}

118

119

sub random_f { # Arguments: ($n,$dfn,$dfd)

120

croak "Usage: random_f(\$n,\$dfn,\$dfd)" if scalar(@_) < 3;

121

my($n, $dfn, $dfd) = @_;

122

croak("($dfn = \$dfn <= 0) or ($dfd = \$dfd <= 0)\nin ".

123

"random_f(\$n,\$dfn,\$dfd)") if (($dfn <= 0) or ($dfd <= 0));

124

return genf($dfn,$dfd) unless wantarray();

125

my $val;

126

my @ans = (0) x $n;

127

foreach $val (@ans) { $val = genf($dfn,$dfd); }

128

return @ans;

129

}

130

131

sub random_gamma { # Arguments: ($n,$a,$r)

132

croak "Usage: random_gamma(\$n,\$a,\$r)" if scalar(@_) < 3;

133

my($n, $a, $r) = @_;

134

croak "($a = \$a <= 0) or ($r = \$r <= 0)\nin random_gamma(\$n,\$a,\$r)"

135

if (($a <= 0) or ($r <= 0));

136

return gengam($a,$r) unless wantarray();

137

my $val;

138

my @ans = (0) x $n;

139

foreach $val (@ans) { $val = gengam($a,$r); }

140

return @ans;

141

}

142

143

sub random_multivariate_normal { # Arguments: ($n, @mean, @covar(2-dim'l))

144

145

croak "Usage: random_multivariate_normal(\$n,\@mean,\@covar(2-dim'l))"

146

if (scalar(@_)) < 3;

147

my $n = shift(@_); # first element is number of obs. desired

148

my $p = scalar(@_)/2; # best guess at dimension of deviate

149

150

# check outline of arguments

151

croak("Sizes of \@mean and \@covar don't match\nin ".

152

"random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")

153

unless (($p == int($p)) and ("$_[$p - 1]" !~ /^ARRAY/) and

154

("$_[$p]" =~ /^ARRAY/));

155

156

# linearize input - it seems faster to push

157

my @linear = ();

158

159

push @linear, splice(@_, 0, $p); # fill first $p slots w/ mean

160

161

# expand array references

162

my $ref;

163

foreach $ref (@_) { # for the rest of the input

164

165

# check length of row of @covariance

166

croak("\@covar is not a $p by $p array ($p is size of \@mean)\nin ".

167

"random_multivariate_normal(\$n, \@mean, \@covar(2-dim'l))")

168

unless (scalar(@{$ref}) == $p);

169

170

push @linear, @{$ref};

171

}

172

173

# load float working array with linearized input

174

putflt(@linear) or

175

croak "Unable to allocate memory\nin random_multivariate_normal";

176

177

# initialize parameter array for multivariate normal generator

178

psetmn($p) or

179

croak "Unable to allocate memory\nin random_multivariate_normal";

180

181

unless (wantarray()) {

182

### if called in a scalar context, returns single refernce to obs

183

pgenmn();

184

return [ getflt($p) ];

185

}

186

187

# otherwise return an $n by $p array of obs.

188

my @ans = (0) x $n;

189

foreach $ref (@ans) {

190

pgenmn();

191

$ref = [ getflt($p) ];

192

}

193

return @ans;

194

}

195

196

sub random_multinomial { # Arguments: ($n,@p)

197

my($n, @p) = @_;

198

my $ncat = scalar(@p); # number of categories

199

$n = int($n);

200

croak "$n = \$n < 0\nin random_multinomial(\$n,\@p)" if ($n < 0);

201

croak "$ncat = (length of \@p) < 2\nin random_multinomial(\$n,\@p)"

202

if ($ncat < 2);

203

rspriw($ncat) or croak "Unable to allocate memory\nin random_multinomial";

204

my($i,$sum,$val) = (0,0,0);

205

pop @p;

206

rsprfw(scalar(@p)) or

207

croak "Unable to allocate memory\nin random_multinomial";

208

foreach $val (@p) {

209

croak "$val = (some \$p[i]) < 0 or > 1\nin random_multinomial(\$n,\@p)"

210

if (($val < 0) or ($val > 1));

211

svprfw($i,$val);

212

$i++;

213

$sum += $val;

214

}

215

croak "Sum of \@p > 1\nin random_multinomial(\$n,\@p)" if ($sum > 0.99999);

216

pgnmul($n, $ncat);

217

### get the results

218

$i = 0;

219

foreach $val (@p) {

220

$val = gvpriw($i);

221

$i++;

222

}

223

push @p, gvpriw($i);

224

return @p;

225

}

226

227

sub random_noncentral_chi_square { # Arguments: ($n,$df,$nonc)

228

croak "Usage: random_noncentral_chi_square(\$n,\$df,\$nonc)"

229

if scalar(@_) < 3;

230

my($n, $df, $nonc) = @_;

231

croak("($df = \$df < 1) or ($nonc = \$nonc) < 0\n".

232

"in random_noncentral_chi_square(\$n,\$df,\$nonc)")

233

if (($df < 1) or ($nonc < 0));

234

return gennch($df,$nonc) unless wantarray();

235

my $val;

236

my @ans = (0) x $n;

237

foreach $val (@ans) { $val = gennch($df,$nonc); }

238

return @ans;

239

}

240

241

sub random_noncentral_f { # Arguments: ($n,$dfn,$dfd,$nonc)

242

croak "Usage: random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)"

243

if scalar(@_) < 4;

244

my($n, $dfn, $dfd, $nonc) = @_;

245

croak("($dfn = \$dfn < 1) or ($dfd = \$dfd <= 0) or ($nonc ".

246

"= \$nonc < 0)\nin random_noncentral_f(\$n,\$dfn,\$dfd,\$nonc)")

247

if (($dfn < 1) or ($dfd <= 0) or ($nonc < 0));

248

return gennf($dfn,$dfd,$nonc) unless wantarray();

249

my $val;

250

my @ans = (0) x $n;

251

foreach $val (@ans) { $val = gennf($dfn,$dfd,$nonc); }

252

return @ans;

253

}

254

255

sub random_normal { # Arguments: ($n,$av,$sd), defaults (1,0,1)

256

return wantarray() ? (gennor(0,1)) : gennor(0,1)

257

if scalar(@_) == 0; # default behavior if no arguments

258

my($n, $av, $sd) = @_;

259

$av = 0 unless defined($av); # $av defaults to 0

260

$sd = 1 unless defined($sd); # $sd defaults to 1, even if $av specified

261

croak "$sd = \$sd < 0\nin random_normal([\$n[,\$av[,\$sd]]])" if ($sd < 0);

262

return gennor($av,$sd) unless wantarray();

263

my $val;

264

my @ans = (0) x $n;

265

foreach $val (@ans) { $val = gennor($av,$sd); }

266

return @ans;

267

}

268

269

sub random_permutation { # Argument: (@array) - array to be permuted.

270

my $n = scalar(@_); # number of elements to be permuted

271

return () if $n == 0;

272

rspriw($n) or

273

croak "Unable to allocate memory\nin random_permutation";

274

pgnprm($n);

275

my($val, $i) = (0,0);

276

my @ans = (0) x $n;

277

foreach $val (@ans) {

278

$val = gvpriw($i);

279

$i++;

280

}

281

return @_[@ans];

282

}

283

284

sub random_permuted_index { # Argument: $n = scalar(@array) (for permutation)

285

croak "Usage: random_permuted_index(\$n)" if scalar(@_) < 1;

286

my $n = int(shift(@_)); # number of elements to be permuted

287

croak "$n = \$n < 0 in random_permuted_index(\$n)" if $n < 0;

288

return () if $n == 0;

289

rspriw($n) or

290

croak "Unable to allocate memory\nin random_permuted_index";

291

pgnprm($n);

292

my($val, $i) = (0,0);

293

my @ans = (0) x $n;

294

foreach $val (@ans) {

295

$val = gvpriw($i);

296

$i++;

297

}

298

return @ans;

299

}

300

301

sub random_uniform { # Arguments: ($n,$low,$high), defaults (1,0,1)

302

return wantarray() ? (genunf(0,1)) : genunf(0,1)

303

if scalar(@_) == 0;

304

croak "Usage: random_uniform([\$n,[\$low,\$high]])"

305

if scalar(@_) == 2; # only default is (0,1) for ($low,$high) both undef

306

my($n, $low, $high) = @_;

307

$low = 0 unless defined($low); # default for $low is 0

308

$high = 1 unless defined($high); # default for $high is 1

309

croak("$low = \$low > \$high = $high\nin ".

310

"random_uniform([\$n,[\$low,\$high]])") if ($low > $high);

311

return genunf($low,$high) unless wantarray();

312

my $val;

313

my @ans = (0) x $n;

314

foreach $val (@ans) { $val = genunf($low,$high); }

315

return @ans;

316

}

317

318

sub random_poisson { # Arguments: ($n, $mu)

319

croak "Usage: random_poisson(\$n,\$mu)" if scalar(@_) < 2;

320

my($n, $mu) = @_;

321

croak "$mu = \$mu < 0\nin random_poisson(\$n,\$mu)" if ($mu < 0);

322

return ignpoi($mu) unless wantarray();

323

my $val;

324

my @ans = (0) x $n;

325

foreach $val (@ans) { $val = ignpoi($mu); }

326

return @ans;

327

}

328

329

sub random_uniform_integer { # Arguments: ($n,$low,$high)

330

croak "Usage: random_uniform_integer(\$n,\$low,\$high)" if scalar(@_) < 3;

331

my($n, $low, $high) = @_;

332

$low = int($low);

333

$high = int($high);

334

croak("$low = \$low > \$high = $high\nin ".

335

"random_uniform_integer(\$n,\$low,\$high)") if ($low > $high);

336

my $range = $high - $low;

337

croak("$range = (\$high - \$low) > 2147483561\nin ".

338

"random_uniform_integer(\$n,\$low,\$high)") if ($range > 2147483561);

339

return ($low + ignuin(0,$range)) unless wantarray();

340

my $val;

341

my @ans = (0) x $n;

342

foreach $val (@ans) { $val = $low + ignuin(0,$range); }

343

return @ans;

344

}

345

346

sub random_negative_binomial { # Arguments: ($n,$ne,$p)

347

croak "Usage: random_negative_binomial(\$n,\$ne,\$p)" if scalar(@_) < 3;

348

my($n, $ne, $p) = @_;

349

$ne = int($ne);

350

croak("($ne = \$ne <= 0) or ($p = \$p <= 0 or >= 1)\nin ".

351

"random_negative_binomial(\$n,\$ne,\$p)")

352

if (($ne <= 0) or (($p <= 0) or ($p >= 1)));

353

return ignnbn($ne,$p) unless wantarray();

354

my $val;

355

my @ans = (0) x $n;

356

foreach $val (@ans) { $val = ignnbn($ne,$p); }

357

return @ans;

358

}

359

360

sub random_binomial { # Arguments: ($n,$nt,$p)

361

croak "Usage: random_binomial(\$n,\$nt,\$p)" if scalar(@_) < 3;

362

my($n, $nt, $p) = @_;

363

$nt = int($nt);

364

croak("($nt = \$nt < 0) or ($p = \$p < 0 or > 1)\nin ".

365

"random_binomial(\$n,\$nt,\$p)")

366

if (($nt < 0) or (($p < 0) or ($p > 1)));

367

return ignbin($nt,$p) unless wantarray();

368

my $val;

369

my @ans = (0) x $n;

370

foreach $val (@ans) { $val = ignbin($nt,$p); }

371

return @ans;

372

}

373

374

#####################################################################

375

# SEED HANDLER FUNCTIONS #

376

#####################################################################

377

378

sub random_seed_from_phrase { # Argument $phrase

379

my $phrase = shift(@_);

380

$phrase ||= "";

381

return phrtsd($phrase);

382

}

383

384

sub random_get_seed { # no argument

385

return getsd();

386

}

387

388

sub random_set_seed_from_phrase { # Argument $phrase

389

my $phrase = shift(@_);

390

$phrase ||= "";

391

salfph($phrase);

392

return 1;

393

}

394

395

sub random_set_seed { # Argument @seed

396

my($seed1,$seed2) = @_;

397

croak("Usage: random_set_seed(\@seed)\n\@seed[0,1] must be two integers ".

398

"in the range (1,1) to (2147483562,2147483398)\nand usually comes ".

399

"from a call to random_get_seed() ".

400

"or\nrandom_seed_from_phrase(\$phrase).")

401

unless (((($seed1 == int($seed1)) and ($seed2 == int($seed2))) and

402

(($seed1 > 0) and ($seed2 > 0))) and

403

(($seed1 < 2147483563) and ($seed2 < 2147483399)));

404

setall($seed1,$seed2);

405

return 1;

406

}

407

408

#####################################################################

409

# HELPER ROUTINES #

410

# These use the C work arrays and are not intended for export #

411

# (Currently only used in random_multivariate_normal) #

412

#####################################################################

413

414

sub getflt {

415

my $n = $_[0];

416

my $val;

417

my $i = 0;

418

my @junk = (0) x $n;

419

foreach $val (@junk) {

420

$val = gvprfw($i);

421

$i++;

422

}

423

return @junk;

424

}

425

426

sub putflt {

427

my $n = scalar(@_);

428

rsprfw($n) or return 0;

429

my $val;

430

my $i = 0;

431

foreach $val (@_) { # load up floats

432

svprfw($i,$val);

433

$i++;

434

}

435

return 1;

436

}

437

438

# Autoload methods go after =cut, and are processed by the autosplit program.

439

440

1;

441

442

__END__

443

444

=head1 NAME

445

446

B<Math::Random> - Random Number Generators

447

448

=head1 SYNOPSIS

449

450

=over 4

451

452

=item *

453

454

use Math::Random;

455

456

Exports the following routines by default (see L<"Default Routines">):

457

458

random_set_seed_from_phrase

459

random_get_seed

460

random_seed_from_phrase

461

random_set_seed

462

random_uniform

463

random_uniform_integer

464

random_permutation

465

random_permuted_index

466

random_normal

467

468

In this case the extended routines (see L<"Extended Routines">) can be

469

used by qualifying them explicitly with C<Math::Random::>, for

470

example: C<$stdexp = Math::Random::random_exponential();>

471

472

=item *

473

474

use Math::Random qw(random_beta

475

random_chi_square

476

random_exponential

477

random_f

478

random_gamma

479

random_multivariate_normal

480

random_multinomial

481

random_noncentral_chi_square

482

random_noncentral_f

483

random_normal

484

random_permutation

485

random_permuted_index

486

random_uniform

487

random_poisson

488

random_uniform_integer

489

random_negative_binomial

490

random_binomial

491

random_seed_from_phrase

492

random_get_seed

493

random_set_seed_from_phrase

494

random_set_seed );

495

496

Exports all the routines explicitly. Use a subset of the list for the

497

routines you want.

498

499

=item *

500

501

use Math::Random qw(:all);

502

503

Exports all the routines, as well.

504

505

=back

506

507

=head1 DESCRIPTION

508

509

B<Math::Random> is a B<Perl> port of the B<C> version of B<randlib>,

510

which is a suite of routines for generating random deviates. See

511

L<"RANDLIB"> for more information.

512

513

This port supports all of the distributions from which the B<Fortran>

514

and B<C> versions generate deviates. The major functionalities that

515

are excluded are the multiple generators/splitting facility and

516

antithetic random number generation. These facilities, along with

517

some of the distributions which I<are> included, are probably not of

518

interest except to the very sophisticated user. If there is

519

sufficient interest, the excluded facilities will be included in a

520

future release. The code to perform the excluded facilities is

521

available as B<randlib> in B<Fortran> and B<C> source.

522

523

=head2 Default Routines

524

525

The routines which are exported by default are the only ones that the

526

average Perl programmer is likely to need.

527

528

=over 4

529

530

=item C<random_set_seed_from_phrase($phrase)>

531

532

Sets the seed of the base generator to a value determined by

533

I<$phrase>. If the module is installed with the default option, the

534

value depends on the machine collating sequence. It should, however,

535

be the same for 7-bit ASCII character strings on all ASCII machines.

536

In the original randlib, the value generated for a given I<$phrase>

537

was consistent from implementation to implementation (it did not rely

538

on the machine collating sequence). Check with your Perl

539

administrator to see if the module was installed with the original

540

seed generator.

541

B<Note:> When the Perl processor loads

542

package B<Math::Random> the seed is set to a value based on the

543

current time. The seed changes each time B<Math::Random> generates

544

something random.

545

546

The ability to set the seed is useful for debugging, or for those who

547

like reproducible runs.

548

549

=item C<random_get_seed()>

550

551

Returns an array of length two which contains the two integers

552

constituting the seed (assuming a call in array context). An

553

invocation in a scalar context returns the integer 2, which is

554

probably not useful.

555

556

=item C<random_seed_from_phrase($phrase)>

557

558

Returns an array of length two which contains the two integers

559

constituting the seed (assuming a call in array context). An

560

invocation in a scalar context returns the integer 2, which is

561

probably not useful. The seed generated is the seed used to set the

562

seed in a call to C<random_set_seed_from_phrase>.

563

564

B<Note:> the following two calls (for the same I<$phrase>) are

565

equivalent:

566

567

random_set_seed(random_seed_from_phrase($phrase));

568

569

and

570

571

random_set_seed_from_phrase($phrase);

572

573

=item C<random_set_seed(@seed)>

574

575

Sets the seed of the base generator to the value I<@seed>[0,1].

576

Usually, the argument I<@seed> should be the result of a call to

577

C<random_get_seed> or C<random_seed_from_phrase>. I<@seed>[0,1] must

578

be two integers in the range S<(1, 1)> to S<(2147483562, 2147483398)>,

579

inclusive.

580

581

=item C<random_uniform($n, $low, $high)>

582

583

=item C<random_uniform($n)>

584

585

=item C<random_uniform()>

586

587

When called in an array context, returns an array of I<$n> deviates

588

generated from a I<uniform($low,>S< >I<$high)> distribution. When

589

called in a scalar context, generates and returns only one such

590

deviate as a scalar, regardless of the value of I<$n>.

591

592

Argument restrictions: I<$low> must be less than or equal to I<$high>.

593

594

Defaults are (1, 0, 1). B<Note:> I<$high> must be specified if

595

I<$low> is specified.

596

597

=item C<random_uniform_integer($n, $low, $high)>

598

599

When called in an array context, returns an array of I<$n> integer

600

deviates generated from a I<uniform($low,>S< >I<$high)> distribution

601

on the integers. When called in a scalar context, generates and

602

returns only one such deviate as a scalar, regardless of the value of

603

I<$n>.

604

605

Argument restrictions: I<$low> and I<$high> are first rounded using

606

C<int()>; the resulting I<$low> must be less than or equal to I<$high>,

607

and the resulting range I<($high - $low)> must not be greater than

608

2147483561.

609

610

There are no defaults; all three arguments must be provided.

611

612

=item C<random_permutation(@array)>

613

614

Returns I<@array>, randomly permuted.

615

616

=item C<random_permuted_index($n)>

617

618

Returns an array of array indices, randomly permuted. The indices

619

used are S<(0, ... ,>(I<$n>S< - >1)). This produces the indices used

620

by C<random_permutation> for a given seed, without passing arrays.

621

622

B<Note:> the following are equivalent:

623

624

random_set_seed_from_phrase('jjv');

625

random_permutation(@array);

626

627

and

628

629

random_set_seed_from_phrase('jjv');

630

@array[(random_permuted_index(scalar(@array)))];

631

632

=item C<random_normal($n, $av, $sd)>

633

634

=item C<random_normal($n, $av)>

635

636

=item C<random_normal($n)>

637

638

=item C<random_normal()>

639

640

When called in an array context, returns an array of I<$n> deviates

641

generated from a I<normal($av, $sd^2)> distribution. When called in a

642

scalar context, generates and returns only one such deviate as a

643

scalar, regardless of the value of I<$n>.

644

645

Argument restrictions: I<$sd> must be non-negative.

646

647

Defaults are (1, 0, 1).

648

649

=back

650

651

=head2 Extended Routines

652

653

These routines generate deviates from many other distributions.

654

655

B<Note:> The parameterizations of these deviates are standard (insofar

656

as there I<is> a standard ... ) but particular attention should be

657

paid to the distributions of the I<beta> and I<gamma> deviates (noted

658

in C<random_beta> and C<random_gamma> below).

659

660

=over 4

661

662

=item C<random_beta($n, $aa, $bb)>

663

664

When called in an array context, returns an array of I<$n> deviates

665

generated from the I<beta> distribution with parameters I<$aa> and

666

I<$bb>. The density of the beta is:

667

668

X^(I<$aa> - 1) * (1 - X)^(I<$bb> - 1) / S(I<$aa> , I<$bb>) for 0 < X <

669

1.

670

671

When called in a scalar context, generates and returns only one such

672

deviate as a scalar, regardless of the value of I<$n>.

673

674

Argument restrictions: Both I<$aa> and I<$bb> must not be less than

675

C<1.0E-37>.

676

677

There are no defaults; all three arguments must be provided.

678

679

=item C<random_binomial($n, $nt, $p)>

680

681

When called in an array context, returns an array of I<$n> outcomes

682

generated from the I<binomial> distribution with number of trials

683

I<$nt> and probability of an event in each trial I<$p>. When called

684

in a scalar context, generates and returns only one such outcome as a

685

scalar, regardless of the value of I<$n>.

686

687

Argument restrictions: I<$nt> is rounded using C<int()>; the result

688

must be non-negative. I<$p> must be between 0 and 1 inclusive.

689

690

There are no defaults; both arguments must be provided.

691

692

=item C<random_chi_square($n, $df)>

693

694

When called in an array context, returns an array of I<$n> deviates

695

generated from the I<chi-square> distribution with I<$df> degrees of

696

freedom. When called in a scalar context, generates and returns only

697

one such deviate as a scalar, regardless of the value of I<$n>.

698

699

Argument restrictions: I<$df> must be positive.

700

701

There are no defaults; both arguments must be provided.

702

703

=item C<random_exponential($n, $av)>

704

705

=item C<random_exponential($n)>

706

707

=item C<random_exponential()>

708

709

When called in an array context, returns an array of I<$n> deviates

710

generated from the I<exponential> distribution with mean I<$av>. When

711

called in a scalar context, generates and returns only one such

712

deviate as a scalar, regardless of the value of I<$n>.

713

714

Argument restrictions: I<$av> must be non-negative.

715

716

Defaults are (1, 1).

717

718

=item C<random_f($n, $dfn, $dfd)>

719

720

When called in an array context, returns an array of I<$n> deviates

721

generated from the I<F> (variance ratio) distribution with degrees of

722

freedom I<$dfn> (numerator) and I<$dfd> (denominator). When called in

723

a scalar context, generates and returns only one such deviate as a

724

scalar, regardless of the value of I<$n>.

725

726

Argument restrictions: Both I<$dfn> and I<$dfd> must be positive.

727

728

There are no defaults; all three arguments must be provided.

729

730

=item C<random_gamma($n, $a, $r)>

731

732

When called in an array context, returns an array of I<$n> deviates

733

generated from the I<gamma> distribution with parameters I<$a> and

734

I<$r>. The density of the gamma is:

735

736

(I<$a>**I<$r>) / Gamma(I<$r>) * X**(I<$r> - 1) * Exp(-I<$a>*X)

737

738

When called in a scalar context, generates and returns only one such

739

deviate as a scalar, regardless of the value of I<$n>.

740

741

Argument restrictions: Both I<$a> and I<$r> must be positive.

742

743

There are no defaults; all three arguments must be provided.

744

745

=item C<random_multinomial($n, @p)>

746

747

When called in an array context, returns single observation from the

748

I<multinomial> distribution, with I<$n> events classified into as many

749

categories as the length of I<@p>. The probability of an event being

750

classified into category I is given by the Ith element of I<@p>.

751

The observation is an array with length equal to I<@p>, so when called

752

in a scalar context it returns the length of @p. The sum of the

753

elements of the observation is equal to I<$n>.

754

755

Argument restrictions: I<$n> is rounded with C<int()> before it is

756

used; the result must be non-negative. I<@p> must have length at

757

least 2. All elements of I<@p> except the last must be between 0 and

758

1 inclusive, and sum to no more than 0.99999. B<Note:> The last

759

element of I<@p> is a dummy to indicate the number of categories, and

760

it is adjusted to bring the sum of the elements of I<@p> to 1.

761

762

There are no defaults; both arguments must be provided.

763

764

=item C<random_multivariate_normal($n, @mean, @covar)>

765

766

When called in an array context, returns an array of I<$n> deviates

767

(each deviate being an array reference) generated from the

768

I<multivariate normal> distribution with mean vector I<@mean> and

769

variance-covariance matrix I<@covar>. When called in a scalar

770

context, generates and returns only one such deviate as an array

771

reference, regardless of the value of I<$n>.

772

773

Argument restrictions: If the dimension of the deviate to be generated

774

is I, I<@mean> should be a length I array of real numbers.

775

I<@covar> should be a length I array of references to length I

776

arrays of real numbers (i.e. a I by I matrix). Further,

777

I<@covar> should be a symmetric positive-definite matrix, although the

778

B<Perl> code does not check positive-definiteness, and the underlying

779

B<C> code assumes the matrix is symmetric. Given that the

780

variance-covariance matrix is symmetric, it doesn't matter if the

781

references refer to rows or columns. If a non-positive definite

782

matrix is passed to the function, it will abort with the following

783

message:

784

785

COVM not positive definite in SETGMN

786

787

Also, a non-symmetric I<@covar> may produce deviates without

788

complaint, although they may not be from the expected distribution.

789

For these reasons, you are encouraged to I<verify the arguments

790

passed>.

791

792

The B<Perl> code I<does> check the dimensionality of I<@mean> and

793

I<@covar> for consistency. It does so by checking that the length of

794

the argument vector passed is odd, that what should be the last

795

element of I<@mean> and the first element of I<@covar> look like they

796

are a number followed by an array reference respectively, and that the

797

arrays referred to in I<@covar> are as long as I<@mean>.

798

799

There are no defaults; all three arguments must be provided.

800

801

=item C<random_negative_binomial($n, $ne, $p)>

802

803

When called in an array context, returns an array of I<$n> outcomes

804

generated from the I<negative binomial> distribution with number of

805

events I<$ne> and probability of an event in each trial I<$p>. When

806

called in a scalar context, generates and returns only one such

807

outcome as a scalar, regardless of the value of I<$n>.

808

809

Argument restrictions: I<$ne> is rounded using C<int()>, the result

810

must be positive. I<$p> must be between 0 and 1 exclusive.

811

812

There are no defaults; both arguments must be provided.

813

814

=item C<random_noncentral_chi_square($n, $df, $nonc)>

815

816

When called in an array context, returns an array of I<$n> deviates

817

generated from the I<noncentral chi-square> distribution with I<$df>

818

degrees of freedom and noncentrality parameter I<$nonc>. When called

819

in a scalar context, generates and returns only one such deviate as a

820

scalar, regardless of the value of I<$n>.

821

822

Argument restrictions: I<$df> must be at least 1, I<$nonc> must be

823

non-negative.

824

825

There are no defaults; all three arguments must be provided.

826

827

=item C<random_noncentral_f($n, $dfn, $dfd, $nonc)>

828

829

When called in an array context, returns an array of I<$n> deviates

830

generated from the I<noncentral F> (variance ratio) distribution with

831

degrees of freedom I<$dfn> (numerator) and I<$dfd> (denominator); and

832

noncentrality parameter I<$nonc>. When called in a scalar context,

833

generates and returns only one such deviate as a scalar, regardless of

834

the value of I<$n>.

835

836

Argument restrictions: I<$dfn> must be at least 1, I<$dfd> must be

837

positive, and I<$nonc> must be non-negative.

838

839

There are no defaults; all four arguments must be provided.

840

841

=item C<random_poisson($n, $mu)>

842

843

When called in an array context, returns an array of I<$n> outcomes

844

generated from the I<Poisson> distribution with mean I<$mu>. When

845

called in a scalar context, generates and returns only one such

846

outcome as a scalar, regardless of the value of I<$n>.

847

848

Argument restrictions: I<$mu> must be non-negative.

849

850

There are no defaults; both arguments must be provided.

851

852

=back

853

854

=head1 ERROR HANDLING

855

856

The B<Perl> code should C<croak> if bad arguments are passed or if the

857

underlying B<C> code cannot allocate the necessary memory. The only

858

error which should kill the job without C<croak>ing is a non-positive

859

definite variance-covariance matrix passed to

860

C<random_multivarite_normal> (see L<"Extended Routines">).

861

862

=head1 RANDLIB

863

864

B<randlib> is available in B<Fortran> and B<C> source form, and will

865

soon be available in B<Fortran90> source as well. B<randlib.c> can be

866

obtained from B<statlib>. Send mail whose message is I<'send

867

randlib.c.shar from general'> to:

868

869

statlib@lib.stat.cmu.edu

870

871

B<randlib.c> can also be obtained by anonymous B<ftp> to:

872

873

odin.mdacc.tmc.edu (143.111.62.32)

874

875

where it is available as

876

877

/pub/source/randlib.c-1.3.tar.gz

878

879

For obvious reasons, the original B<randlib> (in B<Fortran>) has been

880

renamed to

881

882

/pub/source/randlib.f-1.3.tar.gz

883

884

on the same machine.

885

886

Our FTP index is on file C<./pub/index>.

887

888

If you have Internet access and a browser you might note the following

889

web site addresses:

890

891

University of Texas M. D. Anderson Cancer Center Home Page:

892

893

http://www.mdanderson.org/

894

895

Department of Biomathematics Home Page:

896

897

http://odin.mdacc.tmc.edu/

898

899

Available software:

900

901

http://biostatistics.mdanderson.org/SoftwareDownload/

902

903

=head1 SUPPORT

904

905

This work was supported in part by grant CA-16672 from the National

906

Cancer Institute. We are grateful to Larry and Pat McNeil of Corpus

907

Cristi for their generous support. Some equipment used in this effort

908

was provided by IBM as part of a cooperative study agreement; we thank

909

them.

910

911

=head1 CODE MANIPULATION

912

913

The B<C> version of B<randlib> was obtained by translating the

914

original B<Fortran> B<randlib> using B<PROMULA.FORTRAN>, and

915

performing some hand crafting of the result.

916

917

Information on B<PROMULA.FORTRAN> can be obtained from:

918

919

PROMULA Development Corporation

920

3620 N. High Street, Suite 301

921

Columbus, Ohio 43214

922

(614) 263-5454

923

924

F<wrapper.c> (now obsolete) was created by using B<SWIG>, and

925

performing some modification of the result. B<SWIG> also produced the

926

skeleton of F<Random.pm>.

927

928

Information on B<SWIG> can be obtained from:

929

930

http://www.swig.org

931

932

=head1 SOURCES

933

934

The following routines, which were written by others and lightly

935

modified for consistency in packaging, are included in B<randlib>.

936

937

=over 4

938

939

=item Bottom Level Routines

940

941

These routines are a transliteration of the B<Pascal> in the reference

942

to B<Fortran>, and thence to B<C>.

943

944

L'Ecuyer, P., and Cote, S. "Implementing a Random Number Package with

945

Splitting Facilities." ACM Transactions on Mathematical Software,

946

17:98-111 (1991).

947

948

=item Exponential

949

950

This code was obtained from Netlib.

951

952

Ahrens, J. H., and Dieter, U. "Computer Methods for Sampling from the

953

Exponential and Normal Distributions." Comm. ACM, 15,10 (Oct. 1972),

954

873-882.

955

956

=item Gamma

957

958

(Case R >= 1.0)

959

960

Ahrens, J. H., and Dieter, U. "Generating Gamma Variates by a Modified

961

Rejection Technique." Comm. ACM, 25,1 (Jan. 1982), 47-54.

962

Algorithm GD

963

964

(Case 0.0 <= R <= 1.0)

965

966

Ahrens, J. H., and Dieter, U. "Computer Methods for Sampling from

967

Gamma, Beta, Poisson and Binomial Distributions." Computing, 12 (1974),

968

223-246. Adaptation of algorithm GS.

969

970

=item Normal

971

972

This code was obtained from netlib.

973

974

Ahrens, J. H., and Dieter, U. "Extensions of Forsythe's Method for

975

Random Sampling from the Normal Distribution." Math. Comput., 27,124

976

(Oct. 1973), 927-937.

977

978

=item Binomial

979

980

This code was kindly sent to Dr. Brown by Dr. Kachitvichyanukul.

981

982

Kachitvichyanukul, V., and Schmeiser, B. W. "Binomial Random Variate

983

Generation." Comm. ACM, 31, 2 (Feb. 1988), 216.

984

985

=item Poisson

986

987

This code was obtained from netlib.

988

989

Ahrens, J. H., and Dieter, U. "Computer Generation of Poisson Deviates

990

from Modified Normal Distributions." ACM Trans. Math. Software, 8, 2

991

(June 1982), 163-179.

992

993

=item Beta

994

995

This code was written by us following the recipe in the following.

996

997

Cheng, R. C. H. "Generating Beta Variables with Nonintegral Shape

998

Parameters." Comm. ACM, 21:317-322 (1978). (Algorithms BB and BC)

999

1000

=item Linpack

1001

1002

Routines C<SPOFA> and C<SDOT> are used to perform the Cholesky

1003

decomposition of the covariance matrix in C<SETGMN> (used for the

1004

generation of multivariate normal deviates).

1005

1006

Dongarra, J. J., Moler, C. B., Bunch, J. R., and Stewart, G. W.

1007

Linpack User's Guide. SIAM Press, Philadelphia. (1979)

1008

1009

=item Multinomial

1010

1011

The algorithm is from page 559 of Devroye, Luc Non-Uniform Random

1012

Variate Generation. New York: Springer-Verlag, 1986.

1013

1014

=item Negative Binomial

1015

1016

The algorithm is from page 480 of Devroye, Luc Non-Uniform Random

1017

Variate Generation. New York: Springer-Verlag, 1986.

1018

1019

=back

1020

1021

=head1 VERSION

1022

1023

This POD documents B<Math::Random> version 0.71.

1024

1025

=head1 AUTHORS

1026

1027

=over 4

1028

1029

=item *

1030

1031

B<Math::Random> (the B<Perl> port of B<Randlib>) was put together by

1032

John Venier and Barry W. Brown with help from B<SWIG>. For version

1033

0.61, Geoffrey Rommel made various cosmetic changes. Version 0.64 uses

1034

plain vanilla XS rather than SWIG.

1035

1036

=item *

1037

1038

B<randlib> was compiled and written by Barry W. Brown, James Lovato,

1039

Kathy Russell, and John Venier.

1040

1041

=item *

1042

1043

Correspondence regarding B<Math::Random> or B<randlib> should be

1044

addressed to John Venier by email to

1045

1046

jvenier@mdanderson.org

1047

1048

=item *

1049

1050

Our address is:

1051

1052

Department of Biomathematics, Box 237

1053

The University of Texas, M.D. Anderson Cancer Center

1054

1515 Holcombe Boulevard

1055

Houston, TX 77030

1056

1057

=item *

1058

1059

Geoffrey Rommel may be reached at grommel [at] cpan [dot] org.

1060

1061

=back

1062

1063

=head1 LEGALITIES

1064

1065

=over 4

1066

1067

=item *

1068

1069

The programs in the B<Perl> code distributed with B<Math::Random> and

1070

in the B<C> code F<helper.c>, as well as the documentation, are

1071

copyright by John Venier and Barry W. Brown for the University of

1072

Texas M. D. Anderson Cancer Center in 1997. They may be distributed

1073

and used under the same conditions as B<Perl>.

1074

1075

=item *

1076

1077

F<randlib.c>, F<com.c>, and F<randlib.h> are from B<randlib> (See

1078

L<"RANDLIB">) and are distributed with the following legalities.

1079

1080

Code that appeared in an ACM publication is subject to their

1081

algorithms policy:

1082

1083

Submittal of an algorithm for publication in one of the ACM

1084

Transactions implies that unrestricted use of the algorithm within a

1085

computer is permissible. General permission to copy and distribute

1086

the algorithm without fee is granted provided that the copies are not

1087

made or distributed for direct commercial advantage. The ACM

1088

copyright notice and the title of the publication and its date appear,

1089

and notice is given that copying is by permission of the Association

1090

for Computing Machinery. To copy otherwise, or to republish, requires

1091

a fee and/or specific permission.

1092

1093

Krogh, F. "Algorithms Policy." ACM Tran. Math. Softw. 13 (1987),

1094

183-186.

1095

1096

Note, however, that only the particular expression of an algorithm can

1097

be copyrighted, not the algorithm per se; see 17 USC 102E<40>bE<41>.

1098

1099

We place the Randlib code that we have written in the public domain.

1100

1101

=item *

1102

1103

B<Math::Randlib> and B<randlib> are distributed with B<NO WARRANTY>.

1104

See L<"NO WARRANTY">.

1105

1106

=back

1107

1108

=head1 NO WARRANTY

1109

1110

WE PROVIDE ABSOLUTELY NO WARRANTY OF ANY KIND EITHER EXPRESS OR

1111

IMPLIED, INCLUDING BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF

1112

MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK

1113

AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD

1114

THIS PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY

1115

SERVICING, REPAIR OR CORRECTION.

1116

1117

IN NO EVENT SHALL THE UNIVERSITY OF TEXAS OR ANY OF ITS COMPONENT

1118

INSTITUTIONS INCLUDING M. D. ANDERSON HOSPITAL BE LIABLE TO YOU FOR

1119

DAMAGES, INCLUDING ANY LOST PROFITS, LOST MONIES, OR OTHER SPECIAL,

1120

INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR

1121

INABILITY TO USE (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA OR

1122

ITS ANALYSIS BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY THIRD

1123

PARTIES FROM) THE PROGRAM.

1124

1125

(Above NO WARRANTY modified from the GNU NO WARRANTY statement.)

1126

1127

=cut