view release on metacpan or  search on metacpan

lib/Math/Random/AcceptReject.pm  view on Meta::CPAN

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<f(x)> as specified by the C<pdf>
parameter to C<new()>.

=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__

=head1 SEE ALSO

L<http://en.wikipedia.org/wiki/Rejection_sampling>

L<Math::Random::MT>, L<Math::Random>, L<Math::Random::OO>,
L<Math::TrulyRandom>

L<Math::Symbolic>

The examples in the F<examples/> subdirectory of this distribution.

=head1 AUTHOR

Steffen Mueller, E<lt>smueller@cpan.orgE<gt>