File Coverage

blib/lib/Math/Random/AcceptReject.pm

Criterion	Covered	Total	%
statement	15	64	23.4
branch	0	28	0.0
condition	0	6	0.0
subroutine	5	8	62.5
pod	2	2	100.0
total	22	108	20.3

package Math::Random::AcceptReject;

use 5.006;
use strict;
use warnings;

our $VERSION = '0.02';

use Carp qw/croak/;
use Math::Symbolic qw/parse_from_string/;

=head1 NAME

Math::Random::AcceptReject - Acceptance-Rejection PDF transformations

=head1 SYNOPSIS

  use Math::Random::AcceptReject;
  my $pdf = Math::Random::AcceptReject->new(
    xmin   => 0, # defaults to 0
    xmax   => 2, # defaults to 1
    ymax   => 2, # no default!
    pdf    => 'x',    # a triangle from x=0 to x=2 with slope 1
    random => 'rand', # use Perl's builtin (default)
  );

  foreach (1..100) {
      my $rnd = $pdf->rand();
      # ...
  }

  # Use Math::Random::MT instead of bultin rand()
  # Same target PDF but as Perl code instead of a Math::Symbolic
  # function!
  use Math::Random::MT;
  my $mt = Math::Random::Mt->new($seed);
  $pdf = Math::Random::AcceptReject->new(
    xmax   => 2,
    ymax   => 2,
    pdf    => sub { $_[0] },
    random => sub { $mt->rand() };
  );

=head1 DESCRIPTION

This module implements acceptance-rejection transformations of uniformly
distributed random numbers to mostly arbitrary probability density functions
(Is).

Note that whereas J. von Neumann's algorithm can transform from arbitrary
source PDFs to arbitrary destination PDFs, this module is currently limited to
uniform C<[0,1]> source PDFs!

=head1 METHODS

=cut

=head2 new

Creates a new random number generator. Takes named arguments.

Mandatory parameters:

  pdf:  The probability density function. This can either be a
        subroutine reference which takes an argument ('x') and
        returns f(x), a Math::Symbolic tree representing f(x) and
        using the variable 'x', or a string which can be parsed
        as such a Math::Symbolic tree.
  ymax: Maximum value of the target PDF f(x) in the x range. This
        max theoretically be safely set to a very large value which
        is much higher than the real maximum of f(x) within
        the range [xmin,xmax]. The efficiency of the algorithm goes
        down with

Optional parameters:

  random: The random number generator. Defaults to using Perl's
          rand() function. May be set to either 'rand' for the
          default or a subroutine reference for custom random
          number generators. Expected to return one or more(!)
          random numbers per call.
  xmin:   Minimum value for x. Defaults to 0.
  xmax:   Maximum value for x. Defaults to 1.

=cut

sub _dor {
    foreach (@_) {
        return $_ if defined $_;
    }
    return();
}

sub new {
    my $proto = shift;
    my $class = ref($proto)||$proto;

    my %args = @_;

    # Argument checking
    if (not defined $args{ymax}) {
        croak("Need 'ymax' parameter.");
    }
    if (not ref $args{random} eq 'CODE') {
        $args{random} = 'rand';
    }
    if (not defined $args{pdf}) {
        croak("Need 'pdf' parameter.");
    }
    if (not ref($args{pdf})) {
        eval {
            $args{pdf} = parse_from_string($args{pdf});
        };
        if ($@ or not defined $args{pdf}) {
            croak(
                "Error parsing string into Math::Symbolic tree."
                . ($@ ? " Error: $@" : "")
            );
        }
    }
    if (ref($args{pdf}) =~ /^Math::Symbolic/) {
        my ($sub, $leftover) = $args{pdf}->to_sub(x => 0);
        die("Compiling Math::Symbolic tree to sub failed!")
          if not ref($sub) eq 'CODE';
        $args{pdf} = $sub;
    }
    if (not ref($args{pdf}) eq 'CODE') {
        croak("'pdf' parameter needs to be a CODE ref, string, or Math::Symbolic tree.");
    }

    my $self = {
        xmin => _dor($args{xmin}, 0),
        xmax => _dor($args{xmax}, 1),
        ymax => $args{ymax},
        random => $args{random},
        pdf => $args{pdf},
        cache => [],
    };

    if ($self->{xmin} >= $self->{xmax}) {
        croak("'xmin' must be smaller than 'xmax'");
    }

    $self->{xdiff} = $self->{xmax} - $self->{xmin};

    bless $self => $class;

    return $self;
}

=head2 rand

Returns the next random number of PDF C as specified by the C
parameter to C.

=cut

sub rand {
    my $self = shift;
    my $rnd = $self->{random};
    my $pdf = $self->{pdf};
    my $cache = $self->{cache};

    my $accept = 0;
    my $u = 0;
    my $f = -1;
    my $x;
    if (ref($rnd) eq 'CODE') {
        while ($u > $f) {
            push @$cache, $rnd->() if not @$cache;
            $x = $self->{xmin} + shift(@$cache) * $self->{xdiff};

            $f = $pdf->($x);

            push @$cache, $rnd->() if not @$cache;
            $u = shift(@$cache) * $self->{ymax};
        }
    }
    else { # rand
        while ($u > $f) {
            $x = $self->{xmin} + rand() * $self->{xdiff};
            $f = $pdf->($x);
            $u = rand() * $self->{ymax};
        }
    }

    return $x;
}

1;
__END__