1use 5.006; 2use strict; 3use warnings; 4 5package Math::Random::OO::Normal; 6# ABSTRACT: Generates random numbers from the normal (Gaussian) distribution 7our $VERSION = '0.22'; # VERSION 8 9# Required modules 10use Carp; 11use Params::Validate ':all'; 12 13# ISA 14use base qw( Class::Accessor::Fast ); 15 16 17{ 18 my $param_spec = { 19 mean => { type => SCALAR }, 20 stdev => { type => SCALAR } 21 }; 22 23 __PACKAGE__->mk_accessors( keys %$param_spec ); 24 #__PACKAGE__->mk_ro_accessors( keys %$param_spec ); 25 26 sub new { 27 my $class = shift; 28 my $self = bless {}, ref($class) ? ref($class) : $class; 29 if ( @_ > 1 ) { 30 $self->mean( $_[0] ); 31 $self->stdev( abs( $_[1] ) ); 32 } 33 elsif ( @_ == 1 ) { 34 $self->mean( $_[0] ); 35 $self->stdev(1); 36 } 37 else { 38 $self->mean(0); 39 $self->stdev(1); 40 } 41 return $self; 42 } 43} 44 45 46sub seed { 47 my $self = shift; 48 srand( $_[0] ); 49} 50 51 52sub next { 53 my ($self) = @_; 54 my $rnd = rand() || 1e-254; # can't have zero for normals 55 return _ltqnorm($rnd) * $self->stdev + $self->mean; 56} 57 58#--------------------------------------------------------------------------# 59# Function for inverse cumulative normal 60# Used with permission 61# http://home.online.no/~pjacklam/notes/invnorm/impl/acklam/perl/ 62# 63# Input checking removed by DAGOLDEN as the input will be prechecked 64#--------------------------------------------------------------------------# 65 66#<<< No perltidy 67sub _ltqnorm { 68 # Lower tail quantile for standard normal distribution function. 69 # 70 # This function returns an approximation of the inverse cumulative 71 # standard normal distribution function. I.e., given P, it returns 72 # an approximation to the X satisfying P = Pr{Z <= X} where Z is a 73 # random variable from the standard normal distribution. 74 # 75 # The algorithm uses a minimax approximation by rational functions 76 # and the result has a relative error whose absolute value is less 77 # than 1.15e-9. 78 # 79 # Author: Peter J. Acklam 80 # Time-stamp: 2000-07-19 18:26:14 81 # E-mail: pjacklam@online.no 82 # WWW URL: http://home.online.no/~pjacklam 83 84 my $p = shift; 85 # DAGOLDEN: arg checking will be done earlier 86# die "input argument must be in (0,1)\n" unless 0 < $p && $p < 1; 87 88 # Coefficients in rational approximations. 89 my @a = (-3.969683028665376e+01, 2.209460984245205e+02, 90 -2.759285104469687e+02, 1.383577518672690e+02, 91 -3.066479806614716e+01, 2.506628277459239e+00); 92 my @b = (-5.447609879822406e+01, 1.615858368580409e+02, 93 -1.556989798598866e+02, 6.680131188771972e+01, 94 -1.328068155288572e+01 ); 95 my @c = (-7.784894002430293e-03, -3.223964580411365e-01, 96 -2.400758277161838e+00, -2.549732539343734e+00, 97 4.374664141464968e+00, 2.938163982698783e+00); 98 my @d = ( 7.784695709041462e-03, 3.224671290700398e-01, 99 2.445134137142996e+00, 3.754408661907416e+00); 100 101 # Define break-points. 102 my $plow = 0.02425; 103 my $phigh = 1 - $plow; 104 105 # Rational approximation for lower region: 106 if ( $p < $plow ) { 107 my $q = sqrt(-2*log($p)); 108 return ((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) / 109 (((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1); 110 } 111 112 # Rational approximation for upper region: 113 if ( $phigh < $p ) { 114 my $q = sqrt(-2*log(1-$p)); 115 return -((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) / 116 (((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1); 117 } 118 119 # Rational approximation for central region: 120 my $q = $p - 0.5; 121 my $r = $q*$q; 122 return ((((($a[0]*$r+$a[1])*$r+$a[2])*$r+$a[3])*$r+$a[4])*$r+$a[5])*$q / 123 ((((($b[0]*$r+$b[1])*$r+$b[2])*$r+$b[3])*$r+$b[4])*$r+1); 124} 125#>>> 126 1271; 128 129__END__ 130 131=pod 132 133=encoding utf-8 134 135=head1 NAME 136 137Math::Random::OO::Normal - Generates random numbers from the normal (Gaussian) distribution 138 139=head1 VERSION 140 141version 0.22 142 143=head1 SYNOPSIS 144 145 use Math::Random::OO::Normal; 146 push @prngs, 147 Math::Random::OO::Normal->new(), # mean 0, stdev 1 148 Math::Random::OO::Normal->new(5), # mean 5, stdev 1 149 Math::Random::OO::Normal->new(1,3); # mean 1, stdev 3 150 $_->seed(42) for @prngs; 151 print( $_->next() . "\n" ) for @prngs; 152 153=head1 DESCRIPTION 154 155This subclass of L<Math::Random::OO> generates random reals from the normal 156probability distribution, also called the Gaussian or bell-curve distribution. 157 158The module generates random normals from the inverse of the cumulative 159normal distribution using an approximation algorithm developed by Peter J. 160Acklam and released into the public domain. This algorithm claims a 161relative error of less than 1.15e-9 over the entire region. 162 163See http://home.online.no/~pjacklam/notes/invnorm/ for details and discussion. 164 165=head1 METHODS 166 167=head2 C<new> 168 169 $prng1 = Math::Random::OO::Normal->new(); 170 $prng2 = Math::Random::OO::Normal->new($mean); 171 $prng3 = Math::Random::OO::Normal->new($mean,$stdev); 172 173C<new> takes up to two optional parameters and returns a new 174C<Math::Random::OO::Normal> object. With no parameters, the object generates 175random numbers from the standard normal distribution (mean zero, standard 176deviation one). With a single parameter, the object generates random numbers 177with mean equal to the parameter and standard deviation of one. With two 178parameters, the object generates random numbers with mean equal to the first 179parameter and standard deviation equal to the second parameter. (Standard 180deviation should be positive, but this module takes the absolute value of the 181parameter just in case.) 182 183=head2 C<seed> 184 185 $rv = $prng->seed( @seeds ); 186 187This method seeds the random number generator. At the moment, only the 188first seed value matters. It should be a positive integer. 189 190=head2 C<next> 191 192 $rnd = $prng->next(); 193 194This method returns the next random number from the random number generator. 195It does not take any parameters. 196 197=head1 AUTHOR 198 199David Golden <dagolden@cpan.org> 200 201=head1 COPYRIGHT AND LICENSE 202 203This software is Copyright (c) 2013 by David Golden. 204 205This is free software, licensed under: 206 207 The Apache License, Version 2.0, January 2004 208 209=cut 210