File Coverage

blib/lib/Math/Random/Zipf.pm

Criterion	Covered	Total	%
statement	46	46	100.0
branch	6	8	75.0
condition	7	18	38.8
subroutine	10	10	100.0
pod	4	7	57.1
total	73	89	82.0

package Math::Random::Zipf;

use warnings;
use strict;
use POSIX;

our $VERSION = '0.11';

sub new {
    my ($class, $N, $exp) = @_;

    my $self = bless {
	N => $N,
	alpha => $exp,
	cdf => [],
	pmf => [],
    }, ref($class) || $class;

    my $sum = 0;
    for my $k (1 .. $N) {
	$sum += ($self->{pmf}[$k-1] = 1/($k ** $exp));
	push(@{$self->{cdf}}, $sum);
    }
    my $mult = 1/$sum;
    $self->{cdf}->[$_] *= $mult,  $self->{pmf}->[$_] *= $mult for (0 .. $N - 1);

    return $self;
}

sub rand {
    my $self = shift;
    my $x = shift || rand();

    my $tab = $self->{cdf};
    my $lower = -1;
    my $upper = @$tab - 1;
    my $try = -1;
    my $last_try = -1;
    # binary search
    while (1) {
	my $try = POSIX::floor(($lower + $upper + 1) / 2);
	last if $last_try == $try;

	if ($tab->[$try] >= $x) {
	    $upper = $try;
	}
	else {
	    $lower = $try-1;
	}
	$last_try = $try;
    }
    return $upper + 1;
}

sub inv_cdf {
    my ($self, $P) = @_;
    return $self->rand($P);
}

sub pmf {
    my $self = shift;
    my $x = shift;

    return 0 if $x < 1 || $x > $self->{N} || floor($x) != $x;
    return $self->{pmf}->[$x - 1];
}

sub cdf {
    my $self = shift;
    my $x = shift;

    return 0 if $x < 1 || $x > $self->{N} || floor($x) != $x;
    return $self->{cdf}->[$x - 1];
}

sub pmf_ref {
    return shift->{pmf};
}

sub cdf_ref {
    return shift->{cdf};
}


__END__