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